From 3085f3c4c9825035180a8efa8085592cbaf57bf6 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 --- 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 1c2a23a419374..cd6bf5ea512a9 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 d631eb4dad56ae80eaacf901e84693f17edf1c7f 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 --- 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 cd6bf5ea512a9..f832c4323039b 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 6c93bc32f59ecae084a6c34cfe67c68c2145c7e6 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. --- 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 edbe790cbf17d1c3117eae93322d9bddb1734636 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. --- 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 f832c4323039b..0342a58853a94 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 8b399e5febcfa7430d2275d7b4f9f93fc557e72e 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. --- 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; }