From ca5774b3864f61578323bed246647602ca4a77e7 Mon Sep 17 00:00:00 2001 From: Danilo Piparo Date: Tue, 14 Jul 2026 11:52:54 +0200 Subject: [PATCH 1/5] [math][nfc] Improve doc, fix typos (cherry picked from commit 2cebd26e51c9b3af7ccc22cf02974c8b90a300fa) --- math/mathcore/inc/TMath.h | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/math/mathcore/inc/TMath.h b/math/mathcore/inc/TMath.h index 134694e3cf88b..4982c14def89e 100644 --- a/math/mathcore/inc/TMath.h +++ b/math/mathcore/inc/TMath.h @@ -612,7 +612,7 @@ inline Double_t TMath::Tan(Double_t x) { return tan(x); } //////////////////////////////////////////////////////////////////////////////// -/// Returns the hyperbolic sine of `x. +/// Returns the hyperbolic sine of `x`. inline Double_t TMath::SinH(Double_t x) { return sinh(x); } @@ -916,7 +916,7 @@ inline Double_t TMath::QuietNaN() { } //////////////////////////////////////////////////////////////////////////////// -/// Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN). +/// Returns a signaling NaN [as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN). inline Double_t TMath::SignalingNaN() { return std::numeric_limits::signaling_NaN(); @@ -934,7 +934,7 @@ inline Double_t TMath::Infinity() { template inline T TMath::Limits::Min() { - return (std::numeric_limits::min)(); //N.B. use this signature to avoid class with macro min() on Windows + return (std::numeric_limits::min)(); //N.B. use this signature to avoid clashes with macro min() on Windows } //////////////////////////////////////////////////////////////////////////////// @@ -942,7 +942,7 @@ inline T TMath::Limits::Min() { template inline T TMath::Limits::Max() { - return (std::numeric_limits::max)(); //N.B. use this signature to avoid class with macro max() on Windows + return (std::numeric_limits::max)(); //N.B. use this signature to avoid clashes with macro max() on Windows } //////////////////////////////////////////////////////////////////////////////// @@ -1516,7 +1516,7 @@ template Double_t TMath::ModeHalfSample(Long64_t n, const T *a, con const size_t N = std::ceil(n * 0.5); const size_t start = jMin; const size_t stop = start + n - N + 1; // +1 since we use < and not <= - // Find sequentally what v_range is smallest by sliding the half-window + // Find sequentially what v_range is smallest by sliding the half-window for (size_t i = start; i < stop; i++) { Double_t range = values[i + N - 1] - values[i]; if (range < min_v_range) { From 66e69017f0f800b741a800c873eca7ab3a854eaa Mon Sep 17 00:00:00 2001 From: Danilo Piparo Date: Tue, 14 Jul 2026 12:15:39 +0200 Subject: [PATCH 2/5] [math] Fix mistake in Laplacian calculation The standard second-order central, backward and forward differences contained a mistake: a factor-of-4 error. See finite different coefficients here: https://en.wikipedia.org/wiki/Finite_difference_coefficient (cherry picked from commit 99c7a968ecfb4bac70f8bae553f1fe5943378863) --- math/mathcore/inc/TMath.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/math/mathcore/inc/TMath.h b/math/mathcore/inc/TMath.h index 4982c14def89e..23c58962caf92 100644 --- a/math/mathcore/inc/TMath.h +++ b/math/mathcore/inc/TMath.h @@ -1077,15 +1077,15 @@ T *TMath::Laplacian(const Long64_t n, T *f, const double h) T *result = new T[n]; // Forward difference - result[0] = (4 * f[2] + 2 * f[0] - 5 * f[1] - f[3]) / (4 * h * h); + result[0] = (4 * f[2] + 2 * f[0] - 5 * f[1] - f[3]) / (h * h); // Central difference while (i < n - 1) { - result[i] = (f[i + 1] + f[i - 1] - 2 * f[i]) / (4 * h * h); + result[i] = (f[i + 1] + f[i - 1] - 2 * f[i]) / (h * h); i++; } // Backward difference - result[i] = (2 * f[i] - 5 * f[i - 1] + 4 * f[i - 2] - f[i - 3]) / (4 * h * h); + result[i] = (2 * f[i] - 5 * f[i - 1] + 4 * f[i - 2] - f[i - 3]) / (h * h); return result; } From 82a4ef1a48466cebd5d29be9e0b35ea39f0c2556 Mon Sep 17 00:00:00 2001 From: Stephan Hageboeck Date: Wed, 15 Jul 2026 16:12:02 +0200 Subject: [PATCH 3/5] [MathCore] Add a roottest for TMath::Gradient and Laplace. Use a polynomial and its known derivatives to check the correctness of the finite difference methods. (cherry picked from commit 985159909ceff04f1492f94ae90027fc98db3add) --- roottest/root/math/CMakeLists.txt | 2 + roottest/root/math/MathCoreTests.cxx | 77 ++++++++++++++++++++++++++++ 2 files changed, 79 insertions(+) create mode 100644 roottest/root/math/MathCoreTests.cxx diff --git a/roottest/root/math/CMakeLists.txt b/roottest/root/math/CMakeLists.txt index 9a26f3f797c44..1314db5ca11aa 100644 --- a/roottest/root/math/CMakeLists.txt +++ b/roottest/root/math/CMakeLists.txt @@ -1 +1,3 @@ ROOTTEST_ADD_TESTDIRS() + +ROOT_ADD_GTEST(MathCoreTests MathCoreTests.cxx LIBRARIES ROOT::MathCore ROOT::Hist) diff --git a/roottest/root/math/MathCoreTests.cxx b/roottest/root/math/MathCoreTests.cxx new file mode 100644 index 0000000000000..19e85e0eab039 --- /dev/null +++ b/roottest/root/math/MathCoreTests.cxx @@ -0,0 +1,77 @@ +#include +#include + +#include + +#include + +TEST(TMath, Gradient_Laplace) +{ + std::array parameters{2, 100., -1., -2., 0.1}; + std::array parameters1st{}; + std::array parameters2nd{}; + for (unsigned int i = 0; i < 4; ++i) { + parameters1st[i] = parameters[i + 1] * (i + 1); + } + for (unsigned int i = 0; i < 3; ++i) { + parameters2nd[i] = parameters1st[i + 1] * (i + 1); + } + + TF1 poly("poly", "pol4", -10, 10); + TF1 poly1st("derivative", "pol3", -10, 10); + TF1 poly2nd("secondDerivative", "pol2", -10, 10); + poly.SetParameters(parameters.data()); + poly1st.SetParameters(parameters1st.data()); + poly2nd.SetParameters(parameters2nd.data()); + + constexpr std::size_t nPoint = 10000; + std::array vx; + std::array vPoly; + for (unsigned int i = 0; i < nPoint; ++i) { + const auto x = -10. + 20. / nPoint * i; + vx[i] = x; + vPoly[i] = poly.Eval(x); + } + + auto grad = TMath::Gradient(nPoint, vPoly.data(), 20. / nPoint); + auto lap = TMath::Laplacian(nPoint, vPoly.data(), 20. / nPoint); + + auto relativeDiff = [](double val, double ref) { + if (ref == 0.) { + return std::fabs(val - ref); + } + return std::fabs(val - ref) / ref; + }; + + // Check forward/backward differences + EXPECT_LT(relativeDiff(grad[0], poly1st.Eval(-10)), 0.001); + EXPECT_LT(relativeDiff(grad[nPoint - 1], poly1st.Eval(vx[nPoint - 1])), 0.001); + EXPECT_LT(relativeDiff(lap[0], poly2nd.Eval(-10)), 0.001); + EXPECT_LT(relativeDiff(lap[nPoint - 1], poly2nd.Eval(vx[nPoint - 1])), 0.001); + + { + double squaredDiff_grad = 0.; + double maxRelDiff_grad = 0.; + double squaredDiff_laplace = 0.; + double maxRelDiff_laplace = 0.; + // The points on the edges are forward/backward differences, so they will diverge more + // Therefore, run the comparison only on the centre + for (unsigned int i = 1; i < nPoint - 1; ++i) { + const auto x = vx[i]; + const double diff = poly1st.Eval(x) - grad[i]; + squaredDiff_grad += diff * diff; + maxRelDiff_grad = std::max(relativeDiff(grad[i], poly1st.Eval(x)), maxRelDiff_grad); + + const double diff2 = poly2nd.Eval(x) - lap[i]; + squaredDiff_laplace += diff2 * diff2; + maxRelDiff_laplace = std::max(relativeDiff(lap[i], poly2nd.Eval(x)), maxRelDiff_laplace); + } + + // Central differences + EXPECT_LT(maxRelDiff_grad, 0.01); + EXPECT_LT(std::sqrt(squaredDiff_grad), 0.01); + + EXPECT_LT(maxRelDiff_laplace, 0.01); + EXPECT_LT(std::sqrt(squaredDiff_laplace), 0.01); + } +} \ No newline at end of file From 60ff645912ff38dcc1ff52cb434ca74084664bf0 Mon Sep 17 00:00:00 2001 From: Stephan Hageboeck Date: Wed, 15 Jul 2026 18:28:09 +0200 Subject: [PATCH 4/5] [math] Improve interface, docs and errors of finite differences methods. (cherry picked from commit f715d64ae420666092729805638d4494b2cd02fd) --- math/mathcore/inc/TMath.h | 36 +++++++++++++++++++++--------------- 1 file changed, 21 insertions(+), 15 deletions(-) diff --git a/math/mathcore/inc/TMath.h b/math/mathcore/inc/TMath.h index 23c58962caf92..6b66deea74564 100644 --- a/math/mathcore/inc/TMath.h +++ b/math/mathcore/inc/TMath.h @@ -446,8 +446,10 @@ struct Limits { template Iterator LocMax(Iterator first, Iterator last); // Derivatives of an array - template T *Gradient(Long64_t n, T *f, double h = 1); - template T *Laplacian(Long64_t n, T *f, double h = 1); + template + T *Gradient(Long64_t n, T const *f, double h = 1); + template + T *Laplacian(Long64_t n, T const *f, double h = 1); // Hashing ULong_t Hash(const void *txt, Int_t ntxt); @@ -1014,19 +1016,21 @@ Iterator TMath::LocMin(Iterator first, Iterator last) { } //////////////////////////////////////////////////////////////////////////////// -/// \brief Calculate the one-dimensional gradient of an array with length n. +/// \brief Calculate the one-dimensional finite gradient of an array with length n. +/// It is assumed that the values in the array are spaced uniformly in "x", which +/// would amount to taking the derivative w.r.t. x with an infinitesimal step size `h`. +/// /// The first value in the returned array is a forward difference, /// the next n-2 values are central differences, and the last is a backward difference. /// -/// \note Function leads to undefined behavior if n does not match the length of f /// \param n the number of points in the array -/// \param f the array of points. -/// \param h the step size. The default step size is 1. +/// \param f the array of points. It is assumed that these are spaced uniformly in "x". +/// \param h the distance between the points in `f`. /// \return an array of size n with the gradient. Returns nullptr if n < 2 or f empty. Ownership is transferred to the /// caller. template -T *TMath::Gradient(const Long64_t n, T *f, const double h) +T *TMath::Gradient(const Long64_t n, T const *f, const double h) { if (!f) { ::Error("TMath::Gradient", "Input parameter f is empty."); @@ -1052,25 +1056,27 @@ T *TMath::Gradient(const Long64_t n, T *f, const double h) } //////////////////////////////////////////////////////////////////////////////// -/// \brief Calculate the Laplacian of an array with length n. -/// The first value in the returned array is a forward difference, -/// the next n-2 values are central differences, and the last is a backward difference. +/// \brief Calculate second-order finite differences of an array with length n at second-order accuracy. +/// It is assumed that the values in the array are spaced uniformly in "x", where +/// the spacing is represented by the argument `h`. +/// +/// The first value in the returned array is a forward difference at 2nd order accuracy, +/// the next n-2 values are central differences, and the last is the backward difference. /// -/// \note Function leads to undefined behavior if n does not match the length of f /// \param n the number of points in the array -/// \param f the array of points. -/// \param h the step size. The default step size is 1. +/// \param f the array of points. It is assumed that these are spaced uniformly in "x". +/// \param h the distance between the points in `f`. /// \return an array of size n with the laplacian. Returns nullptr if n < 4 or f empty. Ownership is transferred to the /// caller. template -T *TMath::Laplacian(const Long64_t n, T *f, const double h) +T *TMath::Laplacian(const Long64_t n, T const *f, const double h) { if (!f) { ::Error("TMath::Laplacian", "Input parameter f is empty."); return nullptr; } else if (n < 4) { - ::Error("TMath::Laplacian", "Input parameter n=%lld is smaller than 4.", n); + ::Error("TMath::Laplacian", "Need at least four elements, got %lld", n); return nullptr; } Long64_t i = 1; From a068d3a9676625467f2252a0ea3766a00414e795 Mon Sep 17 00:00:00 2001 From: Stephan Hageboeck Date: Thu, 16 Jul 2026 11:36:57 +0200 Subject: [PATCH 5/5] [math] Fix TMath::Gradient/Laplace tests. The test failed without any consequences because it was always returning 0. Other "tests" in this file will need to be updated to check their outputs in the future. Furthermore: - The second derivatives are unstable with short, so remove it from the test. - The targets for the second derivatives were wrong by a factor 4; fix. - Due to rounding, finite differences on integers don't exactly hit the target. Allow a delta of 1. (cherry picked from commit f3d03cb7685d8640c416aa406edf7da37ad91c67) --- math/mathcore/test/testTMath.cxx | 80 ++++++++++++++++++-------------- 1 file changed, 44 insertions(+), 36 deletions(-) diff --git a/math/mathcore/test/testTMath.cxx b/math/mathcore/test/testTMath.cxx index acdc93ff7a134..5ddd07132c03a 100644 --- a/math/mathcore/test/testTMath.cxx +++ b/math/mathcore/test/testTMath.cxx @@ -1,14 +1,14 @@ +#include +#include + +#include #include #include #include #include #include -#include -#include - using std::cout, std::endl, std::vector, std::sort; -using namespace TMath; bool showVector = true; @@ -38,35 +38,45 @@ void testNormCross() } template -void testArrayDerivatives() +bool testArrayDerivatives() { - const Long64_t n = 10; - const double h = 0.1; - T sa[n] = {18, 47, 183, 98, 56, 74, 28, 75, 10, 89}; - T *gradient = TMath::Gradient(n, sa, h); - T *laplacian = TMath::Laplacian(n, sa, h); - - const T gradienta[n] = {290, 825, 255, -635, -120, -140, 5, -90, 70, 790}; - const T laplaciana[n] = {10875, 2675, -5525, 1075, 1500, -1600, 2325, -2800, 3600, 10000}; + bool failure = false; + constexpr Long64_t n = 10; + constexpr double h = 0.1; + constexpr std::array sa = {18, 47, 183, 98, 56, 74, 28, 75, 10, 89}; + T *gradient = TMath::Gradient(sa.size(), sa.data(), h); + T *laplacian = TMath::Laplacian(sa.size(), sa.data(), h); + + // see https://en.wikipedia.org/wiki/Finite_difference_coefficient + constexpr std::array gradienta = {290, 825, 255, -635, -120, -140, 5, -90, 70, 790}; + // central differences are e.g. (sa[x-1] - 2* sa[x] + sa[x+1]) / h^2 + constexpr std::array laplaciana = {43500, 10700, -22100, 4300, 6000, -6400, 9300, -11200, 14400, 40000}; // test results - for (Long64_t i = 0; i < n; i++) { - if (gradient[i] != gradienta[i]) + if (gradient[i] != gradienta[i]) { Error("testArrayDerivatives", "For Gradient, different values found at i = %lld", i); - - if (laplacian[i] != laplaciana[i]) - Error("testArrayDerivatives", "For Laplacian, different values found at i = %lld", i); + failure = true; + } + + // We have to allow a bit of rounding errors, because TMath internally computes in floating-point: + if (abs(static_cast(laplacian[i]) - laplaciana[i]) > 1) { + Error("testArrayDerivatives", "For Laplacian, different values found at i = %lld.", i); + std::cerr << laplacian[i] << " " << laplaciana[i] << "\n"; + failure = true; + } } delete [] gradient; delete [] laplacian; + return failure; } template void testArrayFunctions() { + using namespace TMath; const U n = 10; const U k = 3; U index[n]; @@ -106,6 +116,7 @@ void testArrayFunctions() template void testIteratorFunctions() { + using namespace TMath; const Long64_t n = 10; vector index(n); Long64_t is; @@ -148,9 +159,7 @@ void testPoints(T x, T y) T dx[4] = {0, 0, 2, 2}; T dy[4] = {0, 2, 2, 0}; - cout << "Point(" << x << "," << y << ") IsInside?: " - << IsInside( x, y, n, dx, dy) << endl; - + cout << "Point(" << x << "," << y << ") IsInside?: " << TMath::IsInside(x, y, n, dx, dy) << endl; } template @@ -160,7 +169,7 @@ void testPlane() T dp2[3] = {1,0,0}; T dp3[3] = {0,1,0}; T dn[3]; - Normal2Plane(dp1, dp2, dp3, dn); + TMath::Normal2Plane(dp1, dp2, dp3, dn); cout << "Normal: (" << dn[0] << ", " << dn[1] << ", " @@ -178,7 +187,8 @@ void testBreitWignerRelativistic() for (Int_t i=0;i<=nPoints;i++) { Double_t currentX = xMinimum+i*xStepSize; - cout << "BreitWignerRelativistic(" << currentX << "," << median << "," << gamma << ") = " << BreitWignerRelativistic(currentX,median,gamma) << endl; + cout << "BreitWignerRelativistic(" << currentX << "," << median << "," << gamma + << ") = " << TMath::BreitWignerRelativistic(currentX, median, gamma) << endl; } } @@ -244,8 +254,12 @@ void testHalfSampleMode() R__ASSERT(TMath::ModeHalfSample(testdata7_n, testdata7, weightdata7) == -1); } -void testTMath() +int main() { + // TODO: Warning, many tests don't yet return any success or failure information, + // so don't trust that they will be red. + bool failure = false; + cout << "Starting tests on TMath..." << endl; cout << "\nNormCross tests: " << endl; @@ -265,12 +279,11 @@ void testTMath() cout << "\nArray derivative tests: " << endl; - testArrayDerivatives(); - testArrayDerivatives(); - testArrayDerivatives(); - testArrayDerivatives(); - testArrayDerivatives(); - testArrayDerivatives(); + failure |= testArrayDerivatives(); + failure |= testArrayDerivatives(); + failure |= testArrayDerivatives(); + failure |= testArrayDerivatives(); + failure |= testArrayDerivatives(); cout << "\nIterator functions tests: " << endl; @@ -298,11 +311,6 @@ void testTMath() cout << "\nHalfSampleMode tests: " << endl; testHalfSampleMode(); -} - -int main() -{ - testTMath(); - return 0; + return failure ? 1 : 0; }