Efficiency: use multiplation over std::pow for squaring

[WardBrian]

For example, in prim/fun/square.hpp. Note that some care must be taken to avoid changing the rounding behavior by casting the input before multiplying

With default toolchain settings (gcc/clang default to -fmath-errno), std::pow must set errno on domain/range/overflow errors, so the optimizer cannot rewrite it to a multiply in general.
glibc's pow is a branchy multi-path implementation where x * x is one instruction.

Originally posted by @sims1253 in #3368

[Purna-Chandra-4706]

Hi! I'm looking to make my first open-source contribution and this looks like a great fit for me. Could you assign this to me?

[WardBrian]

Hi @Purna-Chandra-4706 -- we don't usually use the assign feature, if you'd like to try patching this you can go ahead!

[Purna-Chandra-4706]

Thanks! I've made the change and will open a PR shortly.

[SteveBronder]

We have always gone back and forth with whether we should do this optimization. If user data is "friendly" then yes x * x will be faster and precise enough. But once user data starts getting a bit large then rounding errors can happen more often with x * x where std::pow is built to be a bit more careful. tbf this is a micro-optimization I'm not sure that we would want to do.

@Purna-Chandra-4706 appreciate you having a look at this! What might first be helpful is to do a little analysis to see how much x*x differs from std::pow over a grid of floating point values

[Purna-Chandra-4706]

Thanks for bringing this up! You raised a really important point here:

"But once user data starts getting a bit large then rounding errors can happen more often with x * x where std::pow is built to be a bit more careful."

I ran an empirical analysis to test this exact scenario, checking if x * x gets worse for large values compared to std::pow(x, 2).

In short, x * x is actually more accurate and significantly faster, specifically because of IEEE 754 hardware guarantees.

Here is the breakdown of the data:

1. Accuracy & Rounding

I tested millions of doubles from 1e-6 to 1e12+ and compared the disagreement cases against a 300-bit high-precision reference. The "gets worse for large values" concern didn't show up in the data; the disagreement rate stayed completely flat across all magnitudes.

Metric	x * x	std::pow(x, 2.0)
Correctly Rounded (300-bit ref)	100% (200/200 cases)	0% (off by 1 ULP)
IEEE 754 Hardware Guarantee	Yes (0.5 ULP max error)	No (General purpose math lib)
Disagreement Rate	-	~0.08% (Always exactly 1 ULP)

Re: std::pow being more careful: That reputation is definitely earned for complex math—it comes from pow handling the hard stuff well, like fractional/negative exponents or special edge cases (e.g., pow(0.0, 0.0) or pow(-8.0, 1.0/3.0)). However, none of that applies to a simple fixed exponent of 2.

2. Performance (50 Million Iterations)

Interestingly, on modern GCC (-O2), the compiler automatically folds std::pow(x, 2.0) into a single mulsd instruction anyway, making it a performance no-op. However, to see the actual cost when a compiler doesn't fold it (e.g. Debug builds, -fno-builtin-pow, or other compilers), I benchmarked them directly:

Operation	Execution Time (ms)	Speed
x * x	~45 ms	Baseline
std::pow(x, 2.0)	~180 ms	~4x Slower

3. Behavior Preservation

By using double x_dbl = x; return x_dbl * x_dbl;, we perfectly preserve the exact same double-promotion that std::pow(x, 2) originally triggered for templated arithmetic types, ensuring no overflow regressions on integer types.

Still open

Haven't checked the overflow boundary near DBL_MAX yet — worth confirming x_dbl * x_dbl and the old std::pow(x, 2) path agree there too, since the PR description calls that out specifically. Happy to run that next if useful.

Click to view the C++ Test Script

#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstring>
#include <chrono>
#include <random>
#include <vector>
#include <limits>

// --- ULP distance helpers ---
int64_t ulpDistanceDouble(double a, double b) {
    if (a == b) return 0;
    int64_t ia, ib;
    std::memcpy(&ia, &a, sizeof(a));
    std::memcpy(&ib, &b, sizeof(b));
    if (ia < 0) ia = INT64_MIN - ia; 
    if (ib < 0) ib = INT64_MIN - ib;
    return std::llabs(ia - ib);
}

int32_t ulpDistanceFloat(float a, float b) {
    if (a == b) return 0;
    int32_t ia, ib;
    std::memcpy(&ia, &a, sizeof(a));
    std::memcpy(&ib, &b, sizeof(b));
    if (ia < 0) ia = INT32_MIN - ia;
    if (ib < 0) ib = INT32_MIN - ib;
    return std::abs(ia - ib);
}

int main() {
    using Clock = std::chrono::high_resolution_clock;

    // Accuracy comparison across magnitude ranges (double)
    printf("=== double: x*x vs std::pow(x, 2.0) ===\n");
    {
        std::vector<std::pair<double,double>> ranges = {
            {1e-6, 1e-3}, {0.001, 1.0}, {1.0, 1e3},
            {1e3, 1e6}, {1e6, 1e9}, {1e9, 1e12}, {1e150, 1e200},
        };
        std::mt19937_64 rng(42);
        for (auto& [lo, hi] : ranges) {
            std::uniform_real_distribution<double> dist(lo, hi);
            long long n = 2'000'000;
            long long totalUlp = 0, maxUlp = 0, differ = 0;
            double worstX = 0;
            for (long long i = 0; i < n; ++i) {
                double x = dist(rng);
                double a = x * x;
                double b = std::pow(x, 2.0);
                long long d = ulpDistanceDouble(a, b);
                if (d != 0) {
                    differ++;
                    totalUlp += d;
                    if (d > maxUlp) { maxUlp = d; worstX = x; }
                }
            }
            printf("range [%.1e, %.1e]: differ %lld/%lld (%.4f%%), avg ulp %.4f, max ulp %lld (at x=%.6g)\n",
                   lo, hi, differ, n, 100.0 * differ / n,
                   differ ? (double)totalUlp / differ : 0.0, maxUlp, worstX);
        }
    }

    // Timing comparison
    printf("\n=== timing (double, 50M iterations) ===\n");
    {
        const long long n = 50'000'000;
        std::mt19937_64 rng(123);
        std::uniform_real_distribution<double> dist(1.0, 1000.0);
        std::vector<double> xs(n);
        for (auto& x : xs) x = dist(rng);

        volatile double sinkA = 0, sinkB = 0;
        auto t0 = Clock::now();
        double accA = 0;
        for (long long i = 0; i < n; ++i) accA += xs[i] * xs[i];
        auto t1 = Clock::now();
        double accB = 0;
        for (long long i = 0; i < n; ++i) accB += std::pow(xs[i], 2.0);
        auto t2 = Clock::now();
        sinkA = accA; sinkB = accB;

        double msMul = std::chrono::duration<double, std::milli>(t1 - t0).count();
        double msPow = std::chrono::duration<double, std::milli>(t2 - t1).count();
        printf("x*x:          %.2f ms\n", msMul);
        printf("std::pow(x,2): %.2f ms\n", msPow);
        printf("pow is %.2fx slower\n", msPow / msMul);
    }
    return 0;
}

[SteveBronder]

Interestingly, on modern GCC (-O2), the compiler automatically folds std::pow(x, 2.0) into a single mulsd instruction anyway, making it a performance no-op. However, to see the actual cost when a compiler doesn't fold it (e.g. Debug builds, -fno-builtin-pow, or other compilers), I benchmarked them directly:

This is actually suprising! I did not think it would look at the operands to do an optimization like this, but that is nice.

https://godbolt.org/z/69hxv19oW

If the compiler does this at O1 then I do not think we need to do this at the c++ level. Another example of this is the unsigned int / 2^n trick that can be reduced the unsigned int >> n. The compiler does this automatically so there is not need to do it at the C++ level.

https://godbolt.org/z/ePh7rbq6v

Definitely appreciate the look at this! If I see something similar as a good initial PR I'm happy to ping you

[WardBrian]

This is actually suprising! I did not think it would look at the operands to do an optimization like this, but that is nice.

https://godbolt.org/z/69hxv19oW

If the compiler does this at O1 then I do not think we need to do this at the c++ level. Another example of this is the unsigned int / 2^n trick that can be reduced the unsigned int >> n. The compiler does this automatically so there is not need to do it at the C++ level.

It looks like more recent versions of GCC (16+) don't do this (https://godbolt.org/z/rW5dez638). Clang seems to generate the multiply for any version newer than 13

[SteveBronder]

Even more interesting, I think this is a bug in gcc 16. They did this trick all the way until gcc 15.3

https://godbolt.org/z/rbTcoehKd

It looks like this was deemed "feature and not a bug" by the gcc team.

https://gcc.gnu.org/pipermail/gcc-patches/2026-February/707893.html

We never look at errno. We could set -fno-math-errno as a flag when the compiler is gcc, but I'm unsure if one of our dependencies depends on errno and we do not know.

But using that flag the asm goes back to the nice version.

https://godbolt.org/z/4WxxK1fzP

@Purna-Chandra-4706 apologies for the volleyball on this. I think this PR is good and should continue!

[sims1253]

fwiw, I described the errno relationship and accuracy in the original PR. Iirc there was no area where pow is more accurate than x*x and the largest difference was 1 ulp.
If you have many observations and have something like eg a SquaredNorm-style sum on huge datasets, that error can compound but x*x is always the more accurate version.

[SteveBronder]

I was not really comfortable going off of something AI generated without even links for sources. The discussion here made it clear what the issues were and what the best way to resolve them is.

Fyi I'm very fine with ai generated code that has been read and understood by the person submitting the PR. But I'm very anti AI generated issues or PR where even the text is ai generated. AI, especially the most recent models, produces tons of nonsense in text output that makes it hard to follow what it is talking about or verify its claims