diff --git a/.github/workflows/tooling.yml b/.github/workflows/tooling.yml
index 24aa85f2e..36aaa878b 100644
--- a/.github/workflows/tooling.yml
+++ b/.github/workflows/tooling.yml
@@ -39,8 +39,10 @@ jobs:
         OMPI_CXX: clang++
         MPICH_CC: clang
         MPICH_CXX: clang++
-        CXXFLAGS: "-Werror -Wno-error=pass-failed -fsanitize=address,undefined -shared-libsan"
+        CXXFLAGS: "-Werror -Wno-error=pass-failed -Wno-error=self-assign-overloaded -fsanitize=address,undefined -shared-libsan"
         LDFLAGS: "-fsanitize=address,undefined -shared-libsan -fno-omit-frame-pointer"
+        # Note: Clang<=16 triggers on pybind11 code a false-positive of -Wself-assign-overloaded
+        #   https://github.com/llvm/llvm-project/issues/42469
       run: |
         cmake -S . -B build                  \
           -DpyAMReX_IPO=OFF                  \
diff --git a/src/CompensatedReal.H b/src/CompensatedReal.H
new file mode 100644
index 000000000..81b0b3a7d
--- /dev/null
+++ b/src/CompensatedReal.H
@@ -0,0 +1,469 @@
+/* Copyright 2022 The Regents of the University of California, through Lawrence
+*            Berkeley National Laboratory (subject to receipt of any required
+ *           approvals from the U.S. Dept. of Energy). All rights reserved.
+ *
+ * This file is part of ImpactX.
+ *
+ * Authors: Axel Huebl, Chad Mitchell
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef IMPACTX_COMPENSATED_REAL_H
+#define IMPACTX_COMPENSATED_REAL_H
+
+#include <AMReX_REAL.H>
+#include <AMReX_GpuQualifiers.H>
+#include <AMReX_Math.H>
+
+
+namespace impactx
+{
+    /** A compensated floating-point type that automatically uses Kahan-Babuska summation.
+     *
+     * This class wraps a floating-point value and automatically applies the second-order
+     * iterative Kahan-Babuska algorithm for all arithmetic operations to maintain high
+     * precision when accumulating many small values.
+     *
+     * The algorithm follows Klein (2006) to provide high-precision arithmetic that avoids
+     * floating-point precision errors when dealing with large numbers of small values.
+     *
+     * References:
+     * - https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements
+     * - Klein (2006). "A generalized Kahan–Babuška-Summation-Algorithm". in
+     *   Computing. 76 (3–4). Springer-Verlag: 279–293. doi:10.1007/s00607-005-0139-x
+     *
+     * @tparam T Floating-point type (e.g., amrex::Real or amrex::ParticleReal)
+     */
+    template<typename T = amrex::Real>
+    class CompensatedReal
+    {
+        T m_value = T(0.0);     /**< Main value */
+        T m_cs = T(0.0);        /**< First-order compensation for lost low-order bits */
+        T m_ccs = T(0.0);       /**< Second-order compensation for further lost bits */
+
+    public:
+        /** Default constructor */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        explicit CompensatedReal () = default;
+
+        /** Constructor from value with proper compensation reset
+         *
+         * @param value Initial value
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        explicit CompensatedReal (T value) : m_value(value) {
+            reset_compensation();
+        }
+
+        /** Copy constructor
+         *
+         * @param other Other CompensatedReal to copy from
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        explicit CompensatedReal (CompensatedReal const& other) = default;
+
+        /** Assignment operator
+         *
+         * @param other Other CompensatedReal to assign from
+         * @return Reference to this object
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        CompensatedReal& operator=(CompensatedReal const& other) = default;
+
+        /** Assignment operator from value
+         *
+         * @param value Value to assign
+         * @return Reference to this object
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        CompensatedReal& assign (T value) {
+            m_value = value;
+            reset_compensation();
+            return *this;
+        }
+        // note: intentionally commented out to avoid usage like
+        //   s = s + ds
+        // CompensatedReal& operator=(T value) {
+
+        /** Implicit conversion to the underlying type
+         *
+         * @return The compensated value
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        operator T() const {
+            return value();
+        }
+
+        /** Get the final compensated value
+         *
+         * @return The compensated value (main + first-order + second-order compensation)
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        T value () const {
+            return m_value + m_cs + m_ccs;
+        }
+
+        /** Reset compensation variables to zero
+         *
+         * Call this when starting a new tracking sequence to ensure
+         * compensation variables are properly initialized.
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        void reset_compensation () {
+            m_cs = T(0.0);
+            m_ccs = T(0.0);
+        }
+
+        /** Compensated addition operator
+         *
+         * Applies the second-order iterative Kahan-Babuška algorithm
+         * to add a value with high precision.
+         * Automatically promotes other floating-point types to T.
+         *
+         * @param rhs Value to add (any floating-point type)
+         * @return Reference to this object
+         */
+        template<typename U>
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        CompensatedReal& operator+=(U rhs) {
+            T const other = static_cast<T>(rhs);
+
+            // First-order Kahan-Babuška step
+            T const t_sum = m_value + other;
+            T c;
+            if (amrex::Math::abs(m_value) >= amrex::Math::abs(other)) {
+                c = (m_value - t_sum) + other;
+            } else {
+                c = (other - t_sum) + m_value;
+            }
+            m_value = t_sum;
+
+            // Second-order compensation step
+            T const t_comp = m_cs + c;
+            T cc;
+            if (amrex::Math::abs(m_cs) >= amrex::Math::abs(c)) {
+                cc = (m_cs - t_comp) + c;
+            } else {
+                cc = (c - t_comp) + m_cs;
+            }
+            m_cs = t_comp;
+            m_ccs += cc;
+
+            return *this;
+        }
+
+        /** Compensated subtraction operator
+         *
+         * Applies compensated subtraction by adding the negative value.
+         * Automatically promotes other floating-point types to T.
+         *
+         * @param rhs Value to subtract (any floating-point type)
+         * @return Reference to this object
+         */
+        template<typename U>
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        CompensatedReal& operator-=(U rhs) {
+            return operator+=(-rhs);
+        }
+
+        /** Compensated addition with another CompensatedReal
+         *
+         * @param rhs Other CompensatedReal to add
+         * @return Reference to this object
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        CompensatedReal& operator+=(CompensatedReal const& rhs) {
+            // Add the compensated value
+            T const rhs_value = rhs.value();
+            return operator+=(rhs_value);
+        }
+
+        /** Compensated subtraction with another CompensatedReal
+         *
+         * @param rhs Other CompensatedReal to subtract
+         * @return Reference to this object
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        CompensatedReal& operator-=(CompensatedReal const& rhs) {
+            T const rhs_value = rhs.value();
+            return operator+=(-rhs_value);
+        }
+
+        /** Addition operator
+         *
+         * @param rhs Value to add (any floating-point type)
+         * @return New CompensatedReal with result
+         */
+        template<typename U>
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        CompensatedReal operator+(U rhs) const {
+            CompensatedReal result(*this);
+            result += rhs;
+            return CompensatedReal(result);
+        }
+
+        /** Subtraction operator
+         *
+         * @param rhs Value to subtract (any floating-point type)
+         * @return New CompensatedReal with result
+         */
+        template<typename U>
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        CompensatedReal operator-(U rhs) const {
+            CompensatedReal result(*this);
+            result -= rhs;
+            return CompensatedReal(result);
+        }
+
+        /** Addition operator with CompensatedReal
+         *
+         * @param rhs Other CompensatedReal to add
+         * @return New CompensatedReal with result
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        CompensatedReal operator+(CompensatedReal const& rhs) const {
+            CompensatedReal result(*this);
+            result += rhs;
+            return CompensatedReal(result);
+        }
+
+        /** Subtraction operator with CompensatedReal
+         *
+         * @param rhs Other CompensatedReal to subtract
+         * @return New CompensatedReal with result
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        CompensatedReal operator-(CompensatedReal const& rhs) const {
+            CompensatedReal result(*this);
+            result -= rhs;
+            return CompensatedReal(result);
+        }
+
+        /** Equality comparison
+         *
+         * @param rhs Value to compare with
+         * @return True if values are equal
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        bool operator==(T rhs) const {
+            return value() == rhs;
+        }
+
+        /** Inequality comparison
+         *
+         * @param rhs Value to compare with
+         * @return True if values are not equal
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        bool operator!=(T rhs) const {
+            return value() != rhs;
+        }
+
+        /** Less than comparison
+         *
+         * @param rhs Value to compare with
+         * @return True if this value is less than rhs
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        bool operator<(T rhs) const {
+            return value() < rhs;
+        }
+
+        /** Greater than comparison
+         *
+         * @param rhs Value to compare with
+         * @return True if this value is greater than rhs
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        bool operator>(T rhs) const {
+            return value() > rhs;
+        }
+
+        /** Less than or equal comparison
+         *
+         * @param rhs Value to compare with
+         * @return True if this value is less than or equal to rhs
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        bool operator<=(T rhs) const {
+            return value() <= rhs;
+        }
+
+        /** Greater than or equal comparison
+         *
+         * @param rhs Value to compare with
+         * @return True if this value is greater than or equal to rhs
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        bool operator>=(T rhs) const {
+            return value() >= rhs;
+        }
+
+        /** Equality comparison with CompensatedReal
+         *
+         * @param rhs Other CompensatedReal to compare with
+         * @return True if values are equal
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        bool operator==(CompensatedReal const& rhs) const {
+            return value() == rhs.value();
+        }
+
+        /** Inequality comparison with CompensatedReal
+         *
+         * @param rhs Other CompensatedReal to compare with
+         * @return True if values are not equal
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        bool operator!=(CompensatedReal const& rhs) const {
+            return value() != rhs.value();
+        }
+
+        /** Less than comparison with CompensatedReal
+         *
+         * @param rhs Other CompensatedReal to compare with
+         * @return True if this value is less than rhs
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        bool operator<(CompensatedReal const& rhs) const {
+            return value() < rhs.value();
+        }
+
+        /** Greater than comparison with CompensatedReal
+         *
+         * @param rhs Other CompensatedReal to compare with
+         * @return True if this value is greater than rhs
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        bool operator>(CompensatedReal const& rhs) const {
+            return value() > rhs.value();
+        }
+
+        /** Less than or equal comparison with CompensatedReal
+         *
+         * @param rhs Other CompensatedReal to compare with
+         * @return True if this value is less than or equal to rhs
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        bool operator<=(CompensatedReal const& rhs) const {
+            return value() <= rhs.value();
+        }
+
+        /** Greater than or equal comparison with CompensatedReal
+         *
+         * @param rhs Other CompensatedReal to compare with
+         * @return True if this value is greater than or equal to rhs
+         */
+        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+        bool operator>=(CompensatedReal const& rhs) const {
+            return value() >= rhs.value();
+        }
+    };
+
+    /** Convenience typedef for typical usage with particle real type */
+    using CompensatedParticleReal = CompensatedReal<amrex::ParticleReal>;
+
+    // Free function operators for operand flexibility
+    // These allow operations like: 2.5 + cr (where cr is CompensatedReal)
+    // and ensure the result is always CompensatedReal, not the underlying type
+
+    /** Free function addition operator
+     *
+     * @param lhs Left-hand side value
+     * @param rhs Right-hand side CompensatedReal
+     * @return CompensatedReal result
+     */
+    template<typename T>
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    CompensatedReal<T> operator+(T lhs, CompensatedReal<T> const& rhs) {
+        return rhs + lhs;  // Delegate to member operator
+    }
+
+    /** Free function subtraction operator
+     *
+     * @param lhs Left-hand side value
+     * @param rhs Right-hand side CompensatedReal
+     * @return CompensatedReal result
+     */
+    template<typename T>
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    CompensatedReal<T> operator-(T lhs, CompensatedReal<T> const& rhs) {
+        CompensatedReal<T> result(lhs);  // Create CompensatedReal from lhs
+        result -= rhs;                   // Subtract rhs
+        return CompensatedReal<T>(result);
+    }
+
+    /** Free function equality comparison
+     *
+     * @param lhs Left-hand side value
+     * @param rhs Right-hand side CompensatedReal
+     * @return True if values are equal
+     */
+    template<typename T>
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    bool operator==(T lhs, CompensatedReal<T> const& rhs) {
+        return rhs == lhs;  // Delegate to member operator
+    }
+
+    /** Free function inequality comparison
+     *
+     * @param lhs Left-hand side value
+     * @param rhs Right-hand side CompensatedReal
+     * @return True if values are not equal
+     */
+    template<typename T>
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    bool operator!=(T lhs, CompensatedReal<T> const& rhs) {
+        return rhs != lhs;
+    }
+
+    /** Free function less-than comparison
+     *
+     * @param lhs Left-hand side value
+     * @param rhs Right-hand side CompensatedReal
+     * @return True if lhs < rhs
+     */
+    template<typename T>
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    bool operator<(T lhs, CompensatedReal<T> const& rhs) {
+        return rhs > lhs;  // Note: flipped comparison
+    }
+
+    /** Free function greater-than comparison
+     *
+     * @param lhs Left-hand side value
+     * @param rhs Right-hand side CompensatedReal
+     * @return True if lhs > rhs
+     */
+    template<typename T>
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    bool operator>(T lhs, CompensatedReal<T> const& rhs) {
+        return rhs < lhs;
+    }
+
+    /** Free function less-than-or-equal comparison
+     *
+     * @param lhs Left-hand side value
+     * @param rhs Right-hand side CompensatedReal
+     * @return True if lhs <= rhs
+     */
+    template<typename T>
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    bool operator<=(T lhs, CompensatedReal<T> const& rhs) {
+        return rhs >= lhs;
+    }
+
+    /** Free function greater-than-or-equal comparison
+     *
+     * @param lhs Left-hand side value
+     * @param rhs Right-hand side CompensatedReal
+     * @return True if lhs >= rhs
+     */
+    template<typename T>
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    bool operator>=(T lhs, CompensatedReal<T> const& rhs) {
+        return rhs <= lhs;
+    }
+
+} // namespace impactx
+
+#endif // IMPACTX_COMPENSATED_REAL_H
diff --git a/src/elements/CFbend.H b/src/elements/CFbend.H
index a75a4158a..3f33de0e9 100644
--- a/src/elements/CFbend.H
+++ b/src/elements/CFbend.H
@@ -229,7 +229,6 @@ namespace impactx::elements
             amrex::ParticleReal const pz = refpart.pz;
             amrex::ParticleReal const t = refpart.t;
             amrex::ParticleReal const pt = refpart.pt;
-            amrex::ParticleReal const s = refpart.s;
 
             // length of the current slice
             amrex::ParticleReal const slice_ds = m_ds / nslice();
@@ -252,8 +251,8 @@ namespace impactx::elements
             refpart.z = z - (refpart.px - px)/B;
             refpart.t = t - (theta/B)*pt;
 
-            // advance integrated path length
-            refpart.s = s + slice_ds;
+            // advance integrated path length using compensated summation
+            refpart.s += slice_ds;
         }
 
         /** This pushes the covariance matrix. */
diff --git a/src/elements/ChrDrift.H b/src/elements/ChrDrift.H
index 956824609..5ebb8c45f 100644
--- a/src/elements/ChrDrift.H
+++ b/src/elements/ChrDrift.H
@@ -191,7 +191,6 @@ namespace impactx::elements
             amrex::ParticleReal const pz = refpart.pz;
             amrex::ParticleReal const t = refpart.t;
             amrex::ParticleReal const pt = refpart.pt;
-            amrex::ParticleReal const s = refpart.s;
 
             // length of the current slice
             amrex::ParticleReal const slice_ds = m_ds / nslice();
@@ -205,8 +204,8 @@ namespace impactx::elements
             refpart.z = z + step*pz;
             refpart.t = t - step*pt;
 
-            // advance integrated path length
-            refpart.s = s + slice_ds;
+            // advance integrated path length using compensated summation
+            refpart.s += slice_ds;
         }
 
         /** This pushes the covariance matrix. */
diff --git a/src/elements/ChrPlasmaLens.H b/src/elements/ChrPlasmaLens.H
index b50723051..6ee354dda 100644
--- a/src/elements/ChrPlasmaLens.H
+++ b/src/elements/ChrPlasmaLens.H
@@ -276,7 +276,6 @@ namespace impactx::elements
             amrex::ParticleReal const pz = refpart.pz;
             amrex::ParticleReal const t = refpart.t;
             amrex::ParticleReal const pt = refpart.pt;
-            amrex::ParticleReal const s = refpart.s;
 
             // length of the current slice
             amrex::ParticleReal const slice_ds = m_ds / nslice();
@@ -290,8 +289,8 @@ namespace impactx::elements
             refpart.z = z + step*pz;
             refpart.t = t - step*pt;
 
-            // advance integrated path length
-            refpart.s = s + slice_ds;
+            // advance integrated path length using compensated summation
+            refpart.s += slice_ds;
         }
 
         /** This pushes the covariance matrix. */
diff --git a/src/elements/ChrQuad.H b/src/elements/ChrQuad.H
index bc8f69af5..37a855ca9 100644
--- a/src/elements/ChrQuad.H
+++ b/src/elements/ChrQuad.H
@@ -270,7 +270,6 @@ namespace impactx::elements
             amrex::ParticleReal const pz = refpart.pz;
             amrex::ParticleReal const t = refpart.t;
             amrex::ParticleReal const pt = refpart.pt;
-            amrex::ParticleReal const s = refpart.s;
 
             // length of the current slice
             amrex::ParticleReal const slice_ds = m_ds / nslice();
@@ -284,8 +283,8 @@ namespace impactx::elements
             refpart.z = z + step*pz;
             refpart.t = t - step*pt;
 
-            // advance integrated path length
-            refpart.s = s + slice_ds;
+            // advance integrated path length using compensated summation
+            refpart.s += slice_ds;
         }
 
         /** This pushes the covariance matrix. */
diff --git a/src/elements/ChrUniformAcc.H b/src/elements/ChrUniformAcc.H
index 736db4246..81fc710f7 100644
--- a/src/elements/ChrUniformAcc.H
+++ b/src/elements/ChrUniformAcc.H
@@ -249,7 +249,6 @@ namespace impactx::elements
             amrex::ParticleReal const pz = refpart.pz;
             amrex::ParticleReal const t = refpart.t;
             amrex::ParticleReal const pt = refpart.pt;
-            amrex::ParticleReal const s = refpart.s;
 
             // length of the current slice
             amrex::ParticleReal const slice_ds = m_ds / nslice();
@@ -277,8 +276,8 @@ namespace impactx::elements
             refpart.py = py*bgf/bgi;
             refpart.pz = pz*bgf/bgi;
 
-            // advance integrated path length
-            refpart.s = s + slice_ds;
+            // advance integrated path length using compensated summation
+            refpart.s += slice_ds;
         }
 
         /** This pushes the covariance matrix. */
diff --git a/src/elements/ConstF.H b/src/elements/ConstF.H
index 3b33a806d..c068bec7d 100644
--- a/src/elements/ConstF.H
+++ b/src/elements/ConstF.H
@@ -202,7 +202,6 @@ namespace impactx::elements
             amrex::ParticleReal const pz = refpart.pz;
             amrex::ParticleReal const t = refpart.t;
             amrex::ParticleReal const pt = refpart.pt;
-            amrex::ParticleReal const s = refpart.s;
 
             // length of the current slice
             amrex::ParticleReal const slice_ds = m_ds / nslice();
@@ -216,8 +215,8 @@ namespace impactx::elements
             refpart.z = z + step*pz;
             refpart.t = t - step*pt;
 
-            // advance integrated path length
-            refpart.s = s + slice_ds;
+            // advance integrated path length using compensated summation
+            refpart.s += slice_ds;
 
         }
 
diff --git a/src/elements/Drift.H b/src/elements/Drift.H
index 2eeec2268..f5e1d5c12 100644
--- a/src/elements/Drift.H
+++ b/src/elements/Drift.H
@@ -179,7 +179,6 @@ namespace impactx::elements
             amrex::ParticleReal const pz = refpart.pz;
             amrex::ParticleReal const t = refpart.t;
             amrex::ParticleReal const pt = refpart.pt;
-            amrex::ParticleReal const s = refpart.s;
 
             // length of the current slice
             amrex::ParticleReal const slice_ds = m_ds / nslice();
@@ -193,8 +192,8 @@ namespace impactx::elements
             refpart.z = z + step*pz;
             refpart.t = t - step*pt;
 
-            // advance integrated path length
-            refpart.s = s + slice_ds;
+            // advance integrated path length using compensated summation
+            refpart.s += slice_ds;
         }
 
         /** This pushes the covariance matrix. */
diff --git a/src/elements/ExactCFbend.H b/src/elements/ExactCFbend.H
index 08313dbba..a569d10bf 100644
--- a/src/elements/ExactCFbend.H
+++ b/src/elements/ExactCFbend.H
@@ -476,7 +476,6 @@ namespace impactx::elements
             amrex::ParticleReal const pz = refpart.pz;
             amrex::ParticleReal const t = refpart.t;
             amrex::ParticleReal const pt = refpart.pt;
-            amrex::ParticleReal const s = refpart.s;
             amrex::ParticleReal const brho = refpart.rigidity_Tm();
 
 #if AMREX_DEVICE_COMPILE
@@ -522,8 +521,8 @@ namespace impactx::elements
 
             }
 
-            // advance integrated path length
-            refpart.s = s + slice_ds;
+            // advance integrated path length using compensated summation
+            refpart.s += slice_ds;
         }
 
 
diff --git a/src/elements/ExactDrift.H b/src/elements/ExactDrift.H
index 310f1c4de..bc14dcd8e 100644
--- a/src/elements/ExactDrift.H
+++ b/src/elements/ExactDrift.H
@@ -187,7 +187,6 @@ namespace impactx::elements
             amrex::ParticleReal const pz = refpart.pz;
             amrex::ParticleReal const t = refpart.t;
             amrex::ParticleReal const pt = refpart.pt;
-            amrex::ParticleReal const s = refpart.s;
 
             // length of the current slice
             amrex::ParticleReal const slice_ds = m_ds / nslice();
@@ -201,8 +200,8 @@ namespace impactx::elements
             refpart.z = z + step*pz;
             refpart.t = t - step*pt;
 
-            // advance integrated path length
-            refpart.s = s + slice_ds;
+            // advance integrated path length using compensated summation
+            refpart.s += slice_ds;
         }
 
         /** This pushes the covariance matrix. */
diff --git a/src/elements/ExactMultipole.H b/src/elements/ExactMultipole.H
index 2ea2be5d0..3e0014ebb 100644
--- a/src/elements/ExactMultipole.H
+++ b/src/elements/ExactMultipole.H
@@ -416,7 +416,6 @@ namespace impactx::elements
             amrex::ParticleReal const pz = refpart.pz;
             amrex::ParticleReal const t = refpart.t;
             amrex::ParticleReal const pt = refpart.pt;
-            amrex::ParticleReal const s = refpart.s;
 
             // length of the current slice
             amrex::ParticleReal const slice_ds = m_ds / nslice();
@@ -430,8 +429,8 @@ namespace impactx::elements
             refpart.z = z + step*pz;
             refpart.t = t - step*pt;
 
-            // advance integrated path length
-            refpart.s = s + slice_ds;
+            // advance integrated path length using compensated summation
+            refpart.s += slice_ds;
         }
 
         /** This function returns the linear transport map.
diff --git a/src/elements/ExactQuad.H b/src/elements/ExactQuad.H
index d9a98fab8..d2f558f39 100644
--- a/src/elements/ExactQuad.H
+++ b/src/elements/ExactQuad.H
@@ -342,7 +342,6 @@ namespace impactx::elements
             amrex::ParticleReal const pz = refpart.pz;
             amrex::ParticleReal const t = refpart.t;
             amrex::ParticleReal const pt = refpart.pt;
-            amrex::ParticleReal const s = refpart.s;
 
             // length of the current slice
             amrex::ParticleReal const slice_ds = m_ds / nslice();
@@ -356,8 +355,8 @@ namespace impactx::elements
             refpart.z = z + step*pz;
             refpart.t = t - step*pt;
 
-            // advance integrated path length
-            refpart.s = s + slice_ds;
+            // advance integrated path length using compensated summation
+            refpart.s += slice_ds;
         }
 
         /** This function returns the linear transport map.
diff --git a/src/elements/ExactSbend.H b/src/elements/ExactSbend.H
index 3a581cb5d..43232d7a8 100644
--- a/src/elements/ExactSbend.H
+++ b/src/elements/ExactSbend.H
@@ -264,7 +264,6 @@ namespace impactx::elements
             amrex::ParticleReal const pz = refpart.pz;
             amrex::ParticleReal const t = refpart.t;
             amrex::ParticleReal const pt = refpart.pt;
-            amrex::ParticleReal const s = refpart.s;
 
             // length of the current slice
             amrex::ParticleReal const slice_ds = m_ds / nslice();
@@ -302,8 +301,8 @@ namespace impactx::elements
 
             }
 
-            // advance integrated path length
-            refpart.s = s + slice_ds;
+            // advance integrated path length using compensated summation
+            refpart.s += slice_ds;
         }
 
         /** This pushes the covariance matrix. */
diff --git a/src/elements/LinearMap.H b/src/elements/LinearMap.H
index 9dfe3d358..0c67518c7 100644
--- a/src/elements/LinearMap.H
+++ b/src/elements/LinearMap.H
@@ -140,7 +140,6 @@ namespace impactx::elements
                 amrex::ParticleReal const pz = refpart.pz;
                 amrex::ParticleReal const t = refpart.t;
                 amrex::ParticleReal const pt = refpart.pt;
-                amrex::ParticleReal const s = refpart.s;
 
                 // length of the current slice
                 amrex::ParticleReal const slice_ds = m_ds / nslice();
@@ -154,8 +153,8 @@ namespace impactx::elements
                 refpart.z = z + step*pz;
                 refpart.t = t - step*pt;
 
-                // advance integrated path length
-                refpart.s = s + slice_ds;
+                // advance integrated path length using compensated summation
+                refpart.s += slice_ds;
             }
             // else nothing to do for a zero-length element
         }
diff --git a/src/elements/Quad.H b/src/elements/Quad.H
index 36868a935..cfef0aca9 100644
--- a/src/elements/Quad.H
+++ b/src/elements/Quad.H
@@ -218,7 +218,6 @@ namespace impactx::elements
             amrex::ParticleReal const pz = refpart.pz;
             amrex::ParticleReal const t = refpart.t;
             amrex::ParticleReal const pt = refpart.pt;
-            amrex::ParticleReal const s = refpart.s;
 
             // length of the current slice
             amrex::ParticleReal const slice_ds = m_ds / nslice();
@@ -232,8 +231,8 @@ namespace impactx::elements
             refpart.z = z + step*pz;
             refpart.t = t - step*pt;
 
-            // advance integrated path length
-            refpart.s = s + slice_ds;
+            // advance integrated path length using compensated summation
+            refpart.s += slice_ds;
         }
 
         /** This function returns the linear transport map.
diff --git a/src/elements/RFCavity.H b/src/elements/RFCavity.H
index 162c85937..6aeb435fc 100644
--- a/src/elements/RFCavity.H
+++ b/src/elements/RFCavity.H
@@ -282,7 +282,7 @@ namespace RFCavityData
             amrex::ParticleReal const z = refpart.z;
             amrex::ParticleReal const pz = refpart.pz;
             amrex::ParticleReal const pt = refpart.pt;
-            amrex::ParticleReal const s = refpart.s;
+            amrex::ParticleReal const s = refpart.s.value();
             amrex::ParticleReal const sedge = refpart.sedge;
 
             // initialize linear map (deviation) values
@@ -333,8 +333,8 @@ namespace RFCavityData
                }
             }
 
-            // advance integrated path length
-            refpart.s = s + slice_ds;
+            // advance integrated path length using compensated summation
+            refpart.s += slice_ds;
         }
 
         /** This pushes the covariance matrix. */
diff --git a/src/elements/Sbend.H b/src/elements/Sbend.H
index 37c51ae39..1a246e4d9 100644
--- a/src/elements/Sbend.H
+++ b/src/elements/Sbend.H
@@ -239,7 +239,6 @@ namespace impactx::elements
             amrex::ParticleReal const pz = refpart.pz;
             amrex::ParticleReal const t = refpart.t;
             amrex::ParticleReal const pt = refpart.pt;
-            amrex::ParticleReal const s = refpart.s;
 
             // length of the current slice
             amrex::ParticleReal const slice_ds = m_ds / nslice();
@@ -275,8 +274,8 @@ namespace impactx::elements
 
             }
 
-            // advance integrated path length
-            refpart.s = s + slice_ds;
+            // advance integrated path length using compensated summation
+            refpart.s += slice_ds;
 
         }
 
diff --git a/src/elements/SoftQuad.H b/src/elements/SoftQuad.H
index f002de696..84dee2463 100644
--- a/src/elements/SoftQuad.H
+++ b/src/elements/SoftQuad.H
@@ -331,8 +331,8 @@ namespace SoftQuadrupoleData
             refpart.py = py*bgf/bgi;
             refpart.pz = pz*bgf/bgi;
 
-            // advance integrated path length
-            refpart.s = s + slice_ds;
+            // advance integrated path length using compensated summation
+            refpart.s += slice_ds;
         }
 
         /** This pushes the covariance matrix. */
diff --git a/src/elements/SoftSol.H b/src/elements/SoftSol.H
index 653ecbde5..7966292a2 100644
--- a/src/elements/SoftSol.H
+++ b/src/elements/SoftSol.H
@@ -340,8 +340,8 @@ namespace SoftSolenoidData
             refpart.py = py*bgf/bgi;
             refpart.pz = pz*bgf/bgi;
 
-            // advance integrated path length
-            refpart.s = s + slice_ds;
+            // advance integrated path length using compensated summation
+            refpart.s += slice_ds;
         }
 
         /** This pushes the covariance matrix. */
diff --git a/src/elements/Sol.H b/src/elements/Sol.H
index d788ec643..e5b9255b5 100644
--- a/src/elements/Sol.H
+++ b/src/elements/Sol.H
@@ -216,7 +216,6 @@ namespace impactx::elements
             amrex::ParticleReal const pz = refpart.pz;
             amrex::ParticleReal const t = refpart.t;
             amrex::ParticleReal const pt = refpart.pt;
-            amrex::ParticleReal const s = refpart.s;
 
             // length of the current slice
             amrex::ParticleReal const slice_ds = m_ds / nslice();
@@ -230,8 +229,8 @@ namespace impactx::elements
             refpart.z = z + step*pz;
             refpart.t = t - step*pt;
 
-            // advance integrated path length
-            refpart.s = s + slice_ds;
+            // advance integrated path length using compensated summation
+            refpart.s += slice_ds;
         }
 
         /** This pushes the covariance matrix. */
diff --git a/src/elements/diagnostics/BeamMonitor.cpp b/src/elements/diagnostics/BeamMonitor.cpp
index 02d9a0fc5..3244820dd 100644
--- a/src/elements/diagnostics/BeamMonitor.cpp
+++ b/src/elements/diagnostics/BeamMonitor.cpp
@@ -248,13 +248,13 @@ namespace detail {
         auto d_ui = io::Dataset(dtype_ui, {np});
 
         // openPMD 1.* needs "seconds" here, but we fake it as "s"
-        iteration.setTime(ref_part.s);
+        iteration.setTime(ref_part.s.value());
 
         // reference particle information
         beam.setAttribute( "beta_ref", ref_part.beta() );
         beam.setAttribute( "gamma_ref", ref_part.gamma() );
         beam.setAttribute( "beta_gamma_ref", ref_part.beta_gamma() );
-        beam.setAttribute( "s_ref", ref_part.s );
+        beam.setAttribute( "s_ref", ref_part.s.value() );
         beam.setAttribute( "x_ref", ref_part.x );
         beam.setAttribute( "y_ref", ref_part.y );
         beam.setAttribute( "z_ref", ref_part.z );
diff --git a/src/elements/transformation/Insert.cpp b/src/elements/transformation/Insert.cpp
index 0d8fb1b69..9a75a1088 100644
--- a/src/elements/transformation/Insert.cpp
+++ b/src/elements/transformation/Insert.cpp
@@ -8,6 +8,7 @@
  * License: BSD-3-Clause-LBNL
  */
 #include "elements/transformation/Insert.H"
+#include "CompensatedReal.H"
 
 #include <stdexcept>
 #include <type_traits>
@@ -36,8 +37,9 @@ namespace impactx::elements::transformation
 
         std::list<elements::KnownElements> new_list;
 
-        double s = 0.0;  // in meters   // TODO: if we can avoid a global s, we can avoid wasting significant digits for long lattices
-        double s_next_insert = ds;  // in meters
+        // using compensated arithmetic for precision in long sums along s
+        CompensatedReal<double> s(0.0);  // in meters
+        CompensatedReal<double> s_next_insert(ds);  // in meters
 
         while (!list.empty())
         {
@@ -46,7 +48,7 @@ namespace impactx::elements::transformation
             list.pop_front();
 
             // check where the current element ends
-            double cur_s_out;  // in meters
+            CompensatedReal<double> cur_s_out;  // in meters
             std::visit([&s, &cur_s_out](auto &&cur_element)
             {
                 cur_s_out = s + cur_element.ds();
@@ -55,7 +57,7 @@ namespace impactx::elements::transformation
             // case 1: current element is thick and ends after next insert
             if (s_next_insert < cur_s_out)
             {
-                double const s_rel_insert = s_next_insert - s;
+                CompensatedReal<double> const s_rel_insert = s_next_insert - s;
 
                 // split element and shorten each part
                 elements::KnownElements cur_element_leftover = cur_element_variant;
diff --git a/src/particles/CollectLost.cpp b/src/particles/CollectLost.cpp
index d0553002a..0a0006c42 100644
--- a/src/particles/CollectLost.cpp
+++ b/src/particles/CollectLost.cpp
@@ -89,7 +89,7 @@ namespace impactx
         const int s_runtime_index = dest.GetRealCompIndex("s_lost") - dest.NArrayReal;
 
         RefPart const ref_part = source.GetRefParticle();
-        auto const s_lost = ref_part.s;
+        auto const s_lost = ref_part.s.value();
 
         // have to resize here, not in the constructor because grids have not
         // been built when constructor was called.
diff --git a/src/particles/ReferenceParticle.H b/src/particles/ReferenceParticle.H
index 1728dd27d..36c13aa38 100644
--- a/src/particles/ReferenceParticle.H
+++ b/src/particles/ReferenceParticle.H
@@ -10,6 +10,8 @@
 #ifndef IMPACTX_REFERENCE_PARTICLE_H
 #define IMPACTX_REFERENCE_PARTICLE_H
 
+#include "CompensatedReal.H"
+
 #include <ablastr/constant.H>
 
 #include <AMReX_BLassert.H>
@@ -29,7 +31,7 @@ namespace impactx
      */
     struct RefPart
     {
-        amrex::ParticleReal s = 0.0;  ///< integrated orbit path length, in meters
+        CompensatedParticleReal s{0.0};  ///< integrated (compensated sum) orbit path length, in meters
         amrex::ParticleReal x = 0.0;  ///< horizontal position x, in meters
         amrex::ParticleReal y = 0.0;  ///< vertical position y, in meters
         amrex::ParticleReal z = 0.0;  ///< longitudinal position z, in meters
@@ -66,7 +68,7 @@ namespace impactx
             auto old_mass = mass;
             auto old_charge = charge;
 
-            RefPart zero{};
+            RefPart zero;
             *this = zero;
 
             if (keep_mass) { mass = old_mass; }
diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt
index d74fe1cc5..6901a144e 100644
--- a/src/python/CMakeLists.txt
+++ b/src/python/CMakeLists.txt
@@ -4,6 +4,7 @@
 #
 target_sources(pyImpactX
   PRIVATE
+    CompensatedReal.cpp
     distribution.cpp
     elements.cpp
     ImpactX.cpp
diff --git a/src/python/CompensatedReal.cpp b/src/python/CompensatedReal.cpp
new file mode 100644
index 000000000..931b1c7a2
--- /dev/null
+++ b/src/python/CompensatedReal.cpp
@@ -0,0 +1,106 @@
+/* Copyright 2022-2024 The ImpactX Community
+ *
+ * Authors: Axel Huebl, Chad Mitchell
+ * License: BSD-3-Clause-LBNL
+ */
+#include "pyImpactX.H"
+
+#include <CompensatedReal.H>
+#include <AMReX_REAL.H>
+
+#include <cmath>
+
+namespace py = pybind11;
+using namespace impactx;
+
+
+void init_compensatedreal(py::module& m)
+{
+    // Expose the CompensatedParticleReal typedef
+    py::class_<CompensatedParticleReal>(m, "CompensatedParticleReal")
+        .def(py::init<>(),
+             "Default (zero) constructor for compensated ParticleReal type"
+        )
+        .def(py::init<amrex::ParticleReal>(),
+             py::arg("value"),
+             "Constructor from a ParticleReal value"
+        )
+        .def(py::init<CompensatedParticleReal const&>(),
+             py::arg("other"),
+             "Copy constructor"
+        )
+        .def("assign", &CompensatedParticleReal::assign,
+             py::arg("value"),
+             py::return_value_policy::reference_internal,
+             "Assign a new value and reset compensation"
+        )
+
+        // Arithmetic operators
+        .def(py::self += amrex::ParticleReal())
+        .def(py::self -= amrex::ParticleReal())
+        .def(py::self + amrex::ParticleReal())
+        .def(py::self - amrex::ParticleReal())
+        .def(amrex::ParticleReal() + py::self)
+        .def(amrex::ParticleReal() - py::self)
+        .def(py::self += py::self)
+        .def(py::self -= py::self)  // Clang <17 false-positive: -Wself-assign-overloaded
+        .def(py::self + py::self)
+        .def(py::self - py::self)  // Clang <17 false-positive: -Wself-assign-overloaded
+
+        // Comparison operators
+        .def(py::self == amrex::ParticleReal())
+        .def(py::self != amrex::ParticleReal())
+        .def(py::self < amrex::ParticleReal())
+        .def(py::self <= amrex::ParticleReal())
+        .def(py::self > amrex::ParticleReal())
+        .def(py::self >= amrex::ParticleReal())
+        .def(py::self == py::self)
+        .def(py::self != py::self)
+        .def(py::self < py::self)
+        .def(py::self <= py::self)
+        .def(py::self > py::self)
+        .def(py::self >= py::self)
+
+        // Value access
+        .def_property_readonly("value", &CompensatedParticleReal::value,
+             "Get the final compensated value (main + first-order + second-order compensation)"
+        )
+
+        // conversion to float for compatibility
+        .def("__float__", &CompensatedParticleReal::value,
+             "Convert to Python float"
+        )
+        .def("__abs__", [](CompensatedParticleReal const& self) {
+            return std::abs(self.value());
+        }, "Convert to absolute value")
+
+        // Utility methods
+        .def("reset_compensation", &CompensatedParticleReal::reset_compensation,
+             py::return_value_policy::reference_internal,
+             "Reset compensation variables to zero"
+        )
+
+        // String representation
+        .def("__repr__", [](CompensatedParticleReal const& self) {
+            return "CompensatedParticleReal(" + std::to_string(self.value()) + ")";
+        })
+        .def("__str__", [](CompensatedParticleReal const& self) {
+            return std::to_string(self.value());
+        })
+
+        .def_readonly_static("__doc__",
+             "A compensated floating-point type specialized for ParticleReal values.\n"
+             "\n"
+             "This provides high-precision arithmetic for particle tracking calculations\n"
+             "using the second-order iterative Kahan-Babuska algorithm.\n"
+             "\n"
+             "The algorithm follows Klein (2006) to provide high-precision arithmetic that avoids\n"
+             "floating-point precision errors when dealing with large numbers of small values.\n"
+             "\n"
+             "References:\n"
+             "- https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements\n"
+             "- Klein (2006). \"A generalized Kahan–Babuška-Summation-Algorithm\". in\n"
+             "  Computing. 76 (3–4). Springer-Verlag: 279–293. doi:10.1007/s00607-005-0139-x"
+        )
+    ;
+}
diff --git a/src/python/impactx/__init__.py b/src/python/impactx/__init__.py
index d5b5ae647..a246c98cc 100644
--- a/src/python/impactx/__init__.py
+++ b/src/python/impactx/__init__.py
@@ -18,6 +18,7 @@
 
 # import core bindings to C++
 from . import impactx_pybind
+from .compensated import compensated_cumsum, compensated_sum  # noqa
 from .impactx_pybind import *  # noqa
 from .madx_to_impactx import read_beam  # noqa
 
diff --git a/src/python/impactx/compensated.py b/src/python/impactx/compensated.py
new file mode 100644
index 000000000..0859c3579
--- /dev/null
+++ b/src/python/impactx/compensated.py
@@ -0,0 +1,73 @@
+"""
+Compensated arithmetic utilities for ImpactX.
+
+This module provides high-precision arithmetic functions using the CompensatedParticleReal
+type, which implements the Klein (2006) second-order iterative Kahan-Babuška algorithm.
+
+Copyright 2022-2024 The ImpactX Community
+Authors: Axel Huebl, Chad Mitchell
+License: BSD-3-Clause-LBNL
+"""
+
+
+def compensated_sum(values):
+    """
+    Compute the sum of values using compensated arithmetic.
+
+    This function uses CompensatedParticleReal to maintain high precision
+    when summing many small values, avoiding floating-point precision errors.
+
+    Parameters
+    ----------
+    values : iterable
+        An iterable of numeric values to sum
+
+    Returns
+    -------
+    float
+        The compensated sum of all values
+
+    Examples
+    --------
+    >>> compensated_sum([0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1])
+    1.0
+    """
+    from .impactx_pybind import CompensatedParticleReal
+
+    result = CompensatedParticleReal(0.0)
+    for value in values:
+        result += value
+    return result.value
+
+
+def compensated_cumsum(values):
+    """
+    Compute the cumulative sum of values using compensated arithmetic.
+
+    This function uses CompensatedParticleReal to maintain high precision
+    when computing cumulative sums of many small values, avoiding floating-point
+    precision errors that can accumulate in long sequences.
+
+    Parameters
+    ----------
+    values : iterable
+        An iterable of numeric values to compute cumulative sum for
+
+    Returns
+    -------
+    list
+        A list containing the cumulative sum at each position, starting with 0.0
+
+    Examples
+    --------
+    >>> compensated_cumsum([0.1, 0.2, 0.3])
+    [0.0, 0.1, 0.3, 0.6]
+    """
+    from .impactx_pybind import CompensatedParticleReal
+
+    result = CompensatedParticleReal(0.0)
+    cumulative = [0.0]  # Start with initial value
+    for value in values:
+        result += value
+        cumulative.append(result.value)
+    return cumulative
diff --git a/src/python/impactx/plot/Survey.py b/src/python/impactx/plot/Survey.py
index 234a4cf67..5e5d35b79 100644
--- a/src/python/impactx/plot/Survey.py
+++ b/src/python/impactx/plot/Survey.py
@@ -36,9 +36,9 @@ def plot_survey(
     from math import copysign
 
     import matplotlib.pyplot as plt
-    import numpy as np
     from matplotlib.patches import Rectangle
 
+    from ..compensated import compensated_cumsum
     from .ElementColors import get_element_color
 
     charge_qe = 1.0 if ref is None else ref.charge_qe
@@ -47,10 +47,9 @@ def plot_survey(
 
     element_lengths = [element.ds for element in self]
 
-    # NumPy 2.1+ (i.e. Python 3.10+):
-    # element_s = np.cumulative_sum(element_lengths, include_initial=True)
-    # backport:
-    element_s = np.insert(np.cumsum(element_lengths), 0, 0)
+    # Use accurate cumulative sum to avoid floating-point precision errors
+    # when dealing with many small ds values in long lattices
+    element_s = compensated_cumsum(element_lengths)
 
     ax.hlines(0, 0, element_s[-1], color="black", linestyle="--")
 
diff --git a/src/python/pyImpactX.H b/src/python/pyImpactX.H
index d974338dc..de11ed8f1 100644
--- a/src/python/pyImpactX.H
+++ b/src/python/pyImpactX.H
@@ -12,6 +12,7 @@
 #include <pybind11/pybind11.h>
 #include <pybind11/functional.h>
 // include <pybind11/numpy.h>  // not yet used
+#include <pybind11/operators.h>
 #include <pybind11/stl.h>
 #include <pybind11/stl_bind.h>
 
diff --git a/src/python/pyImpactX.cpp b/src/python/pyImpactX.cpp
index 172b1f3ef..271a3f3aa 100644
--- a/src/python/pyImpactX.cpp
+++ b/src/python/pyImpactX.cpp
@@ -20,6 +20,7 @@ using namespace impactx;
 
 
 // forward declarations of exposed classes
+void init_compensatedreal(py::module&);
 void init_distribution(py::module&);
 void init_elements(py::module&);
 void init_ImpactX(py::module&);
@@ -45,6 +46,7 @@ PYBIND11_MODULE(impactx_pybind, m) {
     )pbdoc";
 
     // note: order from parent to child classes
+    init_compensatedreal(m);
     init_distribution(m);
     init_refparticle(m);
     init_impactxparticlecontainer(m);
diff --git a/tests/python/test_lattice_insert.py b/tests/python/test_lattice_insert.py
index 310d4c2aa..d8ff8a899 100644
--- a/tests/python/test_lattice_insert.py
+++ b/tests/python/test_lattice_insert.py
@@ -36,8 +36,8 @@ def test_element_insert():
         filter(lambda el: el.to_dict()["type"] == "BeamMonitor", fodo)
     )
 
-    assert len(fodo) == 63
-    assert len(inserted_monitors) == 29
+    assert len(fodo) == 64
+    assert len(inserted_monitors) == 30
 
     # clean shutdown
     monitor.finalize()