From b04257c25672d89210b284b406eb1d20218e3aca Mon Sep 17 00:00:00 2001 From: Pasukhin Dmitry Date: Sat, 1 Nov 2025 19:56:04 +0000 Subject: [PATCH] Foundation Classes - Move to constexpr Pascal allocator for PLib::Bin (#777) - Converts `BSplCLib::MaxDegree()` to a constexpr function to enable compile-time evaluation - Replaces the dynamically allocated `BinomAllocator` with a template-based constexpr Pascal's triangle that is computed at compile time --- .../TKMath/BSplCLib/BSplCLib.hxx | 2 +- .../TKMath/BSplCLib/BSplCLib.lxx | 2 +- .../TKMath/GTests/PLib_Test.cxx | 124 +++++++++++++++++- src/FoundationClasses/TKMath/PLib/PLib.cxx | 101 +++++++------- 4 files changed, 165 insertions(+), 64 deletions(-) diff --git a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib.hxx b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib.hxx index d1ce304e4c..931bd74bb6 100644 --- a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib.hxx +++ b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib.hxx @@ -349,7 +349,7 @@ public: const Standard_Real Epsilon = 0.0); //! returns the degree maxima for a BSplineCurve. - static Standard_Integer MaxDegree(); + static inline constexpr Standard_Integer MaxDegree(); //! Perform the Boor algorithm to evaluate a point at //! parameter , with and . diff --git a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib.lxx b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib.lxx index c75f7a817d..20c5b815e0 100644 --- a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib.lxx +++ b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib.lxx @@ -20,7 +20,7 @@ //================================================================================================= -inline Standard_Integer BSplCLib::MaxDegree() +inline constexpr Standard_Integer BSplCLib::MaxDegree() { return 25; } diff --git a/src/FoundationClasses/TKMath/GTests/PLib_Test.cxx b/src/FoundationClasses/TKMath/GTests/PLib_Test.cxx index 4809c85aef..5b12e0aece 100644 --- a/src/FoundationClasses/TKMath/GTests/PLib_Test.cxx +++ b/src/FoundationClasses/TKMath/GTests/PLib_Test.cxx @@ -190,20 +190,38 @@ TEST_F(PLibTest, ZeroWeightsHandling) EXPECT_TRUE(true); // We mainly test that no crash occurs } -// Tests for Binomial coefficient function -TEST_F(PLibTest, BinomialCoefficient) +// Tests for Binomial coefficient function - Basic values +TEST_F(PLibTest, BinomialCoefficient_BasicValues) { - // Test known binomial coefficients - EXPECT_NEAR(PLib::Bin(0, 0), 1.0, Precision::Confusion()); + // Test edge cases + EXPECT_NEAR(PLib::Bin(0, 0), 1.0, Precision::Confusion()) << "C(0,0) should be 1"; + + // Test row n=1: [1, 1] + EXPECT_NEAR(PLib::Bin(1, 0), 1.0, Precision::Confusion()); + EXPECT_NEAR(PLib::Bin(1, 1), 1.0, Precision::Confusion()); + + // Test row n=5: [1, 5, 10, 10, 5, 1] EXPECT_NEAR(PLib::Bin(5, 0), 1.0, Precision::Confusion()); - EXPECT_NEAR(PLib::Bin(5, 5), 1.0, Precision::Confusion()); EXPECT_NEAR(PLib::Bin(5, 1), 5.0, Precision::Confusion()); EXPECT_NEAR(PLib::Bin(5, 2), 10.0, Precision::Confusion()); EXPECT_NEAR(PLib::Bin(5, 3), 10.0, Precision::Confusion()); + EXPECT_NEAR(PLib::Bin(5, 4), 5.0, Precision::Confusion()); + EXPECT_NEAR(PLib::Bin(5, 5), 1.0, Precision::Confusion()); + + // Test row n=10: some specific values + EXPECT_NEAR(PLib::Bin(10, 0), 1.0, Precision::Confusion()); + EXPECT_NEAR(PLib::Bin(10, 1), 10.0, Precision::Confusion()); + EXPECT_NEAR(PLib::Bin(10, 2), 45.0, Precision::Confusion()); EXPECT_NEAR(PLib::Bin(10, 3), 120.0, Precision::Confusion()); + EXPECT_NEAR(PLib::Bin(10, 4), 210.0, Precision::Confusion()); + EXPECT_NEAR(PLib::Bin(10, 5), 252.0, Precision::Confusion()); +} +// Tests for Binomial coefficient - Symmetry property +TEST_F(PLibTest, BinomialCoefficient_Symmetry) +{ // Test symmetry property: C(n,k) = C(n,n-k) - for (Standard_Integer n = 1; n <= 10; n++) + for (Standard_Integer n = 0; n <= 20; n++) { for (Standard_Integer k = 0; k <= n; k++) { @@ -213,6 +231,100 @@ TEST_F(PLibTest, BinomialCoefficient) } } +// Tests for Binomial coefficient - Pascal's triangle recurrence +TEST_F(PLibTest, BinomialCoefficient_Recurrence) +{ + // Test Pascal's triangle recurrence: C(n,k) = C(n-1,k-1) + C(n-1,k) + for (Standard_Integer n = 2; n <= 20; n++) + { + for (Standard_Integer k = 1; k < n; k++) + { + Standard_Real expected = PLib::Bin(n - 1, k - 1) + PLib::Bin(n - 1, k); + Standard_Real actual = PLib::Bin(n, k); + EXPECT_NEAR(actual, expected, Precision::Confusion()) + << "Pascal recurrence failed for C(" << n << "," << k << ")"; + } + } +} + +// Tests for Binomial coefficient - Known large values +TEST_F(PLibTest, BinomialCoefficient_LargeValues) +{ + // Test some larger known values + EXPECT_NEAR(PLib::Bin(20, 10), 184756.0, Precision::Confusion()) << "C(20,10)"; + EXPECT_NEAR(PLib::Bin(15, 7), 6435.0, Precision::Confusion()) << "C(15,7)"; + EXPECT_NEAR(PLib::Bin(25, 5), 53130.0, Precision::Confusion()) << "C(25,5)"; + EXPECT_NEAR(PLib::Bin(25, 12), 5200300.0, Precision::Confusion()) << "C(25,12)"; + + // Test maximum supported degree (25) + EXPECT_NEAR(PLib::Bin(25, 0), 1.0, Precision::Confusion()); + EXPECT_NEAR(PLib::Bin(25, 1), 25.0, Precision::Confusion()); + EXPECT_NEAR(PLib::Bin(25, 25), 1.0, Precision::Confusion()); +} + +// Tests for Binomial coefficient - Edge cases with first and last columns +TEST_F(PLibTest, BinomialCoefficient_EdgeColumns) +{ + // Test first column: C(n,0) = 1 for all n + for (Standard_Integer n = 0; n <= 25; n++) + { + EXPECT_NEAR(PLib::Bin(n, 0), 1.0, Precision::Confusion()) << "C(" << n << ",0) should be 1"; + } + + // Test diagonal: C(n,n) = 1 for all n + for (Standard_Integer n = 0; n <= 25; n++) + { + EXPECT_NEAR(PLib::Bin(n, n), 1.0, Precision::Confusion()) + << "C(" << n << "," << n << ") should be 1"; + } + + // Test second column: C(n,1) = n for all n >= 1 + for (Standard_Integer n = 1; n <= 25; n++) + { + EXPECT_NEAR(PLib::Bin(n, 1), Standard_Real(n), Precision::Confusion()) + << "C(" << n << ",1) should be " << n; + } +} + +// Tests for Binomial coefficient - Complete rows verification +TEST_F(PLibTest, BinomialCoefficient_CompleteRows) +{ + // Verify complete row n=6: [1, 6, 15, 20, 15, 6, 1] + const Standard_Real row6[] = {1.0, 6.0, 15.0, 20.0, 15.0, 6.0, 1.0}; + for (Standard_Integer k = 0; k <= 6; k++) + { + EXPECT_NEAR(PLib::Bin(6, k), row6[k], Precision::Confusion()) + << "C(6," << k << ") verification failed"; + } + + // Verify complete row n=8: [1, 8, 28, 56, 70, 56, 28, 8, 1] + const Standard_Real row8[] = {1.0, 8.0, 28.0, 56.0, 70.0, 56.0, 28.0, 8.0, 1.0}; + for (Standard_Integer k = 0; k <= 8; k++) + { + EXPECT_NEAR(PLib::Bin(8, k), row8[k], Precision::Confusion()) + << "C(8," << k << ") verification failed"; + } +} + +// Tests for Binomial coefficient - Sum property +TEST_F(PLibTest, BinomialCoefficient_SumProperty) +{ + // Test sum property: Sum of C(n,k) for k=0 to n equals 2^n + for (Standard_Integer n = 0; n <= 20; n++) + { + Standard_Real sum = 0.0; + Standard_Real expected = std::pow(2.0, n); + + for (Standard_Integer k = 0; k <= n; k++) + { + sum += PLib::Bin(n, k); + } + + EXPECT_NEAR(sum, expected, Precision::Confusion()) + << "Sum property failed for n=" << n << " (sum should be 2^" << n << "=" << expected << ")"; + } +} + // Tests for polynomial evaluation TEST_F(PLibTest, EvalPolynomial) { diff --git a/src/FoundationClasses/TKMath/PLib/PLib.cxx b/src/FoundationClasses/TKMath/PLib/PLib.cxx index 82a5c096dc..603b0653d6 100644 --- a/src/FoundationClasses/TKMath/PLib/PLib.cxx +++ b/src/FoundationClasses/TKMath/PLib/PLib.cxx @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -197,88 +198,76 @@ void PLib::GetPoles(const TColStd_Array1OfReal& FP, } } -// specialized allocator +// specialized allocator for binomial coefficients using compile-time computation namespace { +//! Compile-time Pascal's triangle allocator for binomial coefficients +//! @tparam MaxDegree Maximum degree N for which C(N,P) can be computed +template class BinomAllocator { public: - //! Main constructor - BinomAllocator(const Standard_Integer theMaxBinom) - : myBinom(NULL), - myMaxBinom(theMaxBinom) + //! Constructor - computes Pascal's triangle at compile time + //! Uses the recurrence relation: C(n,k) = C(n-1,k-1) + C(n-1,k) + constexpr BinomAllocator() + : myBinom{} { - Standard_Integer i, im1, ip1, id2, md2, md3, j, k; - Standard_Integer np1 = myMaxBinom + 1; - myBinom = new Standard_Integer*[np1]; - myBinom[0] = new Standard_Integer[1]; - myBinom[0][0] = 1; - for (i = 1; i < np1; ++i) + // Initialize first row: C(0,0) = 1 + myBinom[0][0] = 1; + + // Build Pascal's triangle row by row + for (Standard_Integer i = 1; i <= MaxDegree; ++i) { - im1 = i - 1; - ip1 = i + 1; - id2 = i >> 1; - md2 = im1 >> 1; - md3 = ip1 >> 1; - k = 0; - myBinom[i] = new Standard_Integer[ip1]; - - for (j = 0; j < id2; ++j) - { - myBinom[i][j] = k + myBinom[im1][j]; - k = myBinom[im1][j]; - } - j = id2; - if (j > md2) - j = im1 - j; - myBinom[i][id2] = k + myBinom[im1][j]; + // First and last elements are always 1 + myBinom[i][0] = 1; + myBinom[i][i] = 1; - for (j = ip1 - md3; j < ip1; j++) + // Use recurrence relation for middle elements + for (Standard_Integer j = 1; j < i; ++j) { - myBinom[i][j] = myBinom[i][i - j]; + myBinom[i][j] = myBinom[i - 1][j - 1] + myBinom[i - 1][j]; } } } - //! Destructor - ~BinomAllocator() - { - // free memory - for (Standard_Integer i = 0; i <= myMaxBinom; ++i) - { - delete[] myBinom[i]; - } - delete[] myBinom; - } - - Standard_Real Value(const Standard_Integer N, const Standard_Integer P) const + //! Returns the binomial coefficient C(N,P) + //! @param N the degree (n in C(n,k)) + //! @param P the parameter (k in C(n,k)) + //! @return the value of C(N,P) + //! @note Caller must ensure N and P are within valid range + constexpr Standard_Integer Value(const Standard_Integer N, const Standard_Integer P) const { - Standard_OutOfRange_Raise_if( - N > myMaxBinom, - "PLib, BinomAllocator: requested degree is greater than maximum supported"); - return Standard_Real(myBinom[N][P]); + return myBinom[N][P]; } private: - BinomAllocator(const BinomAllocator&); - BinomAllocator& operator=(const BinomAllocator&); - -private: - Standard_Integer** myBinom; - Standard_Integer myMaxBinom; + Standard_Integer myBinom[MaxDegree + 1][MaxDegree + 1]; }; -// we do not call BSplCLib here to avoid Cyclic dependency detection by WOK -// static BinomAllocator THE_BINOM (BSplCLib::MaxDegree() + 1); -static BinomAllocator THE_BINOM(25 + 1); +//! Thread-safe lazy initialization of compile-time computed binomial coefficients +//! @tparam MaxDegree Maximum degree supported (default 25) +template +inline const BinomAllocator& GetBinomAllocator() +{ + static constexpr BinomAllocator THE_ALLOCATOR{}; + return THE_ALLOCATOR; +} + } // namespace //================================================================================================= Standard_Real PLib::Bin(const Standard_Integer N, const Standard_Integer P) { - return THE_BINOM.Value(N, P); + const auto& aBinom = GetBinomAllocator(); + + Standard_OutOfRange_Raise_if(N < 0 || N > BSplCLib::MaxDegree(), + "PLib::Bin: degree N is out of supported range [0, 25]"); + Standard_OutOfRange_Raise_if(P < 0 || P > N, + "PLib::Bin: parameter P is out of valid range [0, N]"); + + return Standard_Real(aBinom.Value(N, P)); } //================================================================================================= -- 2.39.5