square() for arithmetic types: multiply instead of calling std::pow(x, 2) - #3368
Closed
sims1253 wants to merge 1 commit into
Closed
square() for arithmetic types: multiply instead of calling std::pow(x, 2)#3368sims1253 wants to merge 1 commit into
sims1253 wants to merge 1 commit into
Conversation
…(x, 2) The arithmetic overload of square() calls std::pow(x, 2) even though its own doc comment says the implementation is 'just x * x'. The pow call cannot be constant-folded as well as a multiply and shows up in hot paths: in a GP regression model (gp_exp_quad_cov distances, n=11), square() accounts for 57 pow calls per gradient and 8.9% of gradient instructions; replacing them measured -9.1% Ir/gradient and -13/-15% per-call wall with bit-identical results (glibc). Widen to double first, then multiply: identical to the previous std::pow(x, 2) (a correctly-rounded square equals the rounded product) including the promoted-to-double semantics for integral arguments, where a raw x * x could overflow, and avoiding double-rounding drift for float arguments. On libms without a correctly-rounded pow results may shift by at most 1 ulp. Same treatment for the two scalar-var squared_distance overloads, which used std::pow(a - b, 2) for the value and (a - b) again inside the callback; the difference is now computed once and squared.
Member
|
Please consult the Stan AI policy for contributions at https://github.com/stan-dev/stan/wiki/AI-Contribution-Policy. Fully agentic PRs are not accepted |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
With default toolchain settings (gcc/clang default to
-fmath-errno),std::powmust seterrnoon domain/range/overflow errors, so the optimizer cannot rewrite it to a multiply in general. glibc'spowis a branchy multi-path implementation (~105 instructions per call measured: 3,473,268 Ir over 33,078 calls) wherex * xis one instruction.On the measured GP regression model,
square()accounts for 57powcalls per gradient and 8.9% of gradient instructions (callgrind): 32,889 of 33,078 executedpowcalls come fromgp_exp_quad_cov(55 kernel pairs for N = 11, plussquare(sigma)andsquare(l)).Widen to
doublefirst, then multiply. The template is enabled for all arithmeticT; a naivex * xwould compute an int product for integralx(overflow for |x| > 46,341 where the promoted path does not) and round twice forfloatx. Widening first reproduces the previous promote-then-round semantics exactly, minus the libm callEval
GP regression model (
gp_exp_quad_cov, n = 11 → 55 kernel pairs), matched binaries, identical inputs, gcc 16.2.1, glibc, Zen 3. Medians of 3 interleaved reps; callgrind on a fixed seeded run (warmup 50 / samples 50, 577 gradient calls in both arms):Checklist
Copyright holder: Maximilian Scholz
The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
- Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
- Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)
the basic tests are passing
./runTests.py test/unit)make test-headers)make test-math-dependencies)make doxygen)make cpplint)the code is written in idiomatic C++ and changes are documented in the doxygen
the new changes are tested