Skip to content

Conversation

@ax3l
Copy link
Member

@ax3l ax3l commented Oct 21, 2025

Long sums, e.g., s of the reference particle, suffer dramatically from floating point precision errors as large and small terms are added in succession thousands to millions of times. Using a compensated sum can keep these numerical errors in check.

This implements, with some vibe coding assistance, the 2nd-order method of Klein (2006) that is linked on Wikipedia. The cost is that s is now of the size of 3 instead of 1 double, e.g., when copying to kernels. We could also use a first order method, which only needs 2 doubles (e.g., Neumaier (1974)).

The latter might be worth considering if we need to use compensated sums for further properties, e.g., x,y,z,t(,px,py,pz,pt) of the reference particle...?

Originally inspired by #1183, but needed at many locations in ImpactX, e.g., existing survey plots #1040, inserting elements #877 #1095 etc. Exposed to Python and used there, too.

To Do

  • define a test, tracking errors in s, x, pt, t, etc. of the reference particle -- .e.g, LHC size and beam lifetime
  • decide if 1st (Neumann) or 2nd order (Klein) is sensible to use

@ax3l ax3l requested a review from cemitch99 October 21, 2025 05:54
@ax3l ax3l force-pushed the s-compensated branch 5 times, most recently from 9204323 to f4e9203 Compare October 21, 2025 06:59

// openPMD 1.* needs "seconds" here, but we fake it as "s"
iteration.setTime(ref_part.s);
iteration.setTime(ref_part.s.value());
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Interesting. Should we just use

Suggested change
iteration.setTime(ref_part.s.value());
iteration.setTime(ref_part.t.value());

?

(Should be a separate PR, just noticed here...)

@ax3l ax3l force-pushed the s-compensated branch 4 times, most recently from ad48adf to cee4e07 Compare October 21, 2025 07:46
*/
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
bool operator==(CompensatedReal const& rhs) const {
return value() == rhs.value();

Check notice

Code scanning / CodeQL

Equality test on floating-point values Note

Equality checks on floating point values can yield unexpected results.
@ax3l ax3l changed the title Compensated Sum for Global Coordinates [WIP] Compensated Sum for Global Coordinates Oct 23, 2025
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants