From cd8d07dc4d7a82c8ef07717a334ad0b1ecd7856b Mon Sep 17 00:00:00 2001 From: silviana amethyst Date: Thu, 16 Jul 2026 07:32:26 +0000 Subject: [PATCH 1/2] fix(pseg): max_cycle_number is a ceiling, not a floor; clamp before unsigned conversion ComputeBoundOnCycleNumber applied the max_cycle_number config with max() instead of min(), so the 'largest cycle number to consider' never bounded anything: a near-unity sample ratio (slow convergence -- high multiplicity, or a slow diverger) made the amplified estimate itself the bound, costing hundreds of Hermite solves per approximation -- and a ratio within ~1e-14 of 1 amplifies past UINT_MAX, where the unclamped conversion to unsigned is undefined behavior. Clamp in the real type before converting: anything not provably below the ceiling (including inf/NaN) becomes the ceiling. The old expectation 'max(5,6) = 6' in compute_bound_on_cycle_num pinned the inverted behavior; it is now 5 (the amplified estimate, under the ceiling). Named regression test covers both the slow-ratio cap and the UB-magnitude clamp. Co-Authored-By: Claude Fable 5 --- .../include/bertini2/endgames/powerseries.hpp | 17 +++-- core/test/endgames/generic_pseg_test.hpp | 72 ++++++++++++++++++- 2 files changed, 82 insertions(+), 7 deletions(-) diff --git a/core/include/bertini2/endgames/powerseries.hpp b/core/include/bertini2/endgames/powerseries.hpp index e5c3e822b..3ed71190d 100644 --- a/core/include/bertini2/endgames/powerseries.hpp +++ b/core/include/bertini2/endgames/powerseries.hpp @@ -378,13 +378,22 @@ class PowerSeriesEndgame : RealT estimate = log(static_cast(this->EndgameSettings().sample_factor))/log(abs(rand_sum2/rand_sum1)); - if (estimate < 1) // would be nan if sample points are same as each other + const auto& ps_config = this->template Get(); + if (estimate < 1) upper_bound_on_cycle_number_ = 1; else { - using std::max; - auto upper_bound = unsigned(round(estimate)*this->template Get().cycle_number_amplification); - upper_bound_on_cycle_number_ = max(upper_bound,this->template Get().max_cycle_number); + // max_cycle_number is a CEILING on the candidate search (each candidate costs a full + // Hermite solve). Clamp before any conversion to unsigned: near-unity sample ratios + // (slow convergence -- high multiplicity, or a slow diverger) drive the estimate toward + // +inf, and unsigned(inf) is undefined behavior. The !(a < b) form also routes NaN to + // the ceiling. This was `max` -- the ceiling was a floor, and the search was unbounded. + using std::round; + RealT amplified = round(estimate) * static_cast(ps_config.cycle_number_amplification); + if (!(amplified < static_cast(ps_config.max_cycle_number))) + upper_bound_on_cycle_number_ = ps_config.max_cycle_number; + else + upper_bound_on_cycle_number_ = std::max(1u, static_cast(amplified)); } return upper_bound_on_cycle_number_; diff --git a/core/test/endgames/generic_pseg_test.hpp b/core/test/endgames/generic_pseg_test.hpp index 62601fd75..b845c3f4e 100644 --- a/core/test/endgames/generic_pseg_test.hpp +++ b/core/test/endgames/generic_pseg_test.hpp @@ -298,7 +298,7 @@ BOOST_AUTO_TEST_CASE(compute_bound_on_cycle_num) my_endgame.ComputeBoundOnCycleNumber(); - BOOST_CHECK(my_endgame.UpperBoundOnCycleNumber() == 6); // max_cycle_num implemented max(5,6) = 6 + BOOST_CHECK(my_endgame.UpperBoundOnCycleNumber() == 5); // round(estimate)*amplification = 5, under the ceiling of 6 [[maybe_unused]] auto first_upper_bound = my_endgame.UpperBoundOnCycleNumber(); @@ -319,9 +319,75 @@ BOOST_AUTO_TEST_CASE(compute_bound_on_cycle_num) my_endgame.ComputeBoundOnCycleNumber(); - BOOST_CHECK(my_endgame.UpperBoundOnCycleNumber() == 6); // max_cycle_num implemented max(5,6) = 6 + BOOST_CHECK(my_endgame.UpperBoundOnCycleNumber() == 5); // round(estimate)*amplification = 5, under the ceiling of 6 -} // end compute bound on cycle number +} // end compute bound on cycle number + + +/** +Regression: max_cycle_number is a CEILING on the cycle-number candidate search. It was applied +with max() instead of min(), so a near-unity sample ratio (slow convergence -- high multiplicity, +or a slow diverger) made the amplified estimate the bound: hundreds of candidates, each costing a +full Hermite solve -- and a ratio close enough to 1 makes the amplified estimate exceed UINT_MAX, +where the old unclamped conversion to unsigned was undefined behavior. +*/ +BOOST_AUTO_TEST_CASE(cycle_number_upper_bound_capped_for_near_unity_sample_ratios) +{ + DefaultPrecision(ambient_precision); + + bertini::System sys; + Var x = Variable::Make("x"); + sys.AddFunction(pow(x-1,3)); + + VariableGroup vars{x}; + sys.AddVariableGroup(vars); + + auto precision_config = PrecisionConfig(sys); + TrackerType tracker(sys); + + bertini::tracking::SteppingConfig stepping_settings; + bertini::tracking::NewtonConfig newton_settings; + tracker.Setup(TestedPredictor, 1e-5, 1e5, stepping_settings, newton_settings); + tracker.PrecisionSetup(precision_config); + + bertini::TimeCont times; + bertini::SampCont samples; + Vec sample(1); + + // consecutive sample differences shrink by a ratio of 0.99: the estimate + // log(sample_factor)/log(0.99) ~ 69, amplified ~ 345 -- far above the ceiling + times.push_back(ComplexFromString(".1")); + sample << ComplexFromString("1.0"); samples.push_back(sample); + times.push_back(ComplexFromString(".05")); + sample << ComplexFromString("2.0"); samples.push_back(sample); // diff 1 + times.push_back(ComplexFromString(".025")); + sample << ComplexFromString("2.99"); samples.push_back(sample); // diff 0.99 + + bertini::endgame::EndgameConfig endgame_settings; + TestedEGType my_endgame(tracker, endgame_settings); + my_endgame.SetTimes(times); + my_endgame.SetSamples(samples); + my_endgame.SetRandVec(1); + + my_endgame.ComputeBoundOnCycleNumber(); + auto ceiling = bertini::endgame::PowerSeriesConfig().max_cycle_number; + BOOST_CHECK_EQUAL(my_endgame.UpperBoundOnCycleNumber(), ceiling); + + // ratio within 1e-14 of 1: the amplified estimate is ~3.5e14 > UINT_MAX -- the old + // unclamped conversion to unsigned was undefined behavior; the bound must be the ceiling + samples.clear(); times.clear(); + times.push_back(ComplexFromString(".1")); + sample << ComplexFromString("1.0"); samples.push_back(sample); + times.push_back(ComplexFromString(".05")); + sample << ComplexFromString("2.0"); samples.push_back(sample); // diff 1 + times.push_back(ComplexFromString(".025")); + sample << ComplexFromString("2.99999999999999"); samples.push_back(sample); // diff 1 - 1e-14 + + my_endgame.SetTimes(times); + my_endgame.SetSamples(samples); + my_endgame.ComputeBoundOnCycleNumber(); + BOOST_CHECK_EQUAL(my_endgame.UpperBoundOnCycleNumber(), ceiling); +} // end cycle_number_upper_bound_capped_for_near_unity_sample_ratios From 1f113eb59b8983fad63c3bbfab0391968b785efc Mon Sep 17 00:00:00 2001 From: silviana amethyst Date: Thu, 16 Jul 2026 07:44:02 +0000 Subject: [PATCH 2/2] fix(pseg): HermiteInterpolateAndSolve now evaluates the actual Hermite interpolant The Horner reconstruction walked the doubled node list at half speed (node z_ii paired with coefficient a_{2*ii}), evaluating a polynomial that was NOT the interpolant: a cubic -- which 3-node Hermite must reproduce exactly -- came back 1.6e-5 off. The limit as the sample window slides to the target was unaffected, so PSEG still converged, but at degraded order: extra iterations per path, poorer approximations than the samples support. The pinned test expectation ('found using matlab') was itself not the Hermite value, and its loose tolerance let the wrong code pass. All expectations re-derived in exact rational arithmetic: the interpolant of the x^8+1 test data at 0 is exactly 0.999999998837890625, and the window-halving errors are 1.162e-9 -> 4.539e-12 -> 1.773e-14 (~256x per halving -- the restored order, visible in-test). New named regression test pins cubic exactness at 1e-20. Computational-semantics note (records doctrine): PSEG approximations change value with this fix -- this is the behavior-epoch case; the epochs mechanism is not yet built, so this note is the marker. Co-Authored-By: Claude Fable 5 --- .../bertini2/endgames/interpolation.hpp | 27 +++++------ core/test/endgames/generic_interpolation.hpp | 10 ++-- core/test/endgames/generic_pseg_test.hpp | 47 +++++++++++++++++-- 3 files changed, 62 insertions(+), 22 deletions(-) diff --git a/core/include/bertini2/endgames/interpolation.hpp b/core/include/bertini2/endgames/interpolation.hpp index 613964461..8cffa2835 100644 --- a/core/include/bertini2/endgames/interpolation.hpp +++ b/core/include/bertini2/endgames/interpolation.hpp @@ -112,20 +112,19 @@ template } } - //Start of Result from Hermite polynomial, this is using the diagonal of the - //finite difference matrix. - Vec Result = space_differences(2*num_sample_points - 1,2*num_sample_points - 1); - - - //This builds the hermite polynomial from the highest term down. - //As we multiply the previous result we will construct the highest term down to the last term. - for (unsigned ii=num_sample_points-1; ii >= 1; --ii) - { - Result = ((Result*(target_time - time_differences(ii)) + space_differences(2*ii, 2*ii)) * (target_time - time_differences(ii-1)) + space_differences(2*ii-1, 2*ii-1)).eval(); - } - - // Last term in hermite polynomial. - return (Result * (target_time - time_differences(0)) + space_differences(0,0)).eval(); + //The interpolant in Newton form is P(x) = sum_k a_k prod_{j Result = space_differences(2*num_sample_points - 1,2*num_sample_points - 1); + for (int k = 2*static_cast(num_sample_points) - 2; k >= 0; --k) + Result = (Result*(target_time - time_differences(k)) + space_differences(k,k)).eval(); + return Result; } //re: HermiteInterpolateAndSolve }} // re: namespaces diff --git a/core/test/endgames/generic_interpolation.hpp b/core/test/endgames/generic_interpolation.hpp index c30da56cc..459df07f9 100644 --- a/core/test/endgames/generic_interpolation.hpp +++ b/core/test/endgames/generic_interpolation.hpp @@ -223,8 +223,12 @@ BOOST_AUTO_TEST_CASE(eight_degree_univariate_advanced_gets_better) Vec< BCT > third_approx = HermiteInterpolateAndSolve(target_time,num_samples,times,samples,derivatives); - BOOST_CHECK((first_approx - correct).norm() < 1e-10); - BOOST_CHECK((second_approx - correct).norm() < 1e-10); - BOOST_CHECK((third_approx - correct).norm() < 1e-10); + // Tolerances calibrated to the TRUE Hermite interpolation errors of these windows (exact + // rational arithmetic): 1.162e-9, 4.539e-12, 1.773e-14 -- shrinking ~256x per halving. + // The old flat 1e-10 was calibrated to a mis-indexed Horner that happened to land closer + // to the truth than the actual interpolant does on the first window. + BOOST_CHECK((first_approx - correct).norm() < 2e-9); + BOOST_CHECK((second_approx - correct).norm() < 1e-11); + BOOST_CHECK((third_approx - correct).norm() < 1e-13); }//end hermite test case \ No newline at end of file diff --git a/core/test/endgames/generic_pseg_test.hpp b/core/test/endgames/generic_pseg_test.hpp index b845c3f4e..0e59b2cfb 100644 --- a/core/test/endgames/generic_pseg_test.hpp +++ b/core/test/endgames/generic_pseg_test.hpp @@ -114,12 +114,45 @@ BOOST_AUTO_TEST_CASE( basic_hermite_test_case_against_matlab ) Vec< BCT > first_approx = HermiteInterpolateAndSolve(target_time,num_samples,times,samples,derivatives); - BOOST_CHECK( norm(first_approx(0) - ComplexFromString("0.9999999767578209232082898114211261253459","0")) < 1e-7); - // answer was found using matlab for a check. difference is diff is 2.32422e-08 + // The unique Hermite interpolant of this data, evaluated at 0, is EXACTLY the terminating + // decimal below (computed independently in exact rational arithmetic). The previous + // expectation (0.99999997675782..., attributed to matlab) was not the Hermite value, and the + // loose 1e-7 tolerance let a mis-indexed Horner reconstruction pass against it. + BOOST_CHECK( norm(first_approx(0) - ComplexFromString("0.999999998837890625","0")) < 1e-20); }//end basic hermite test case mp against matlab +/** +Regression: Hermite interpolation with n nodes (values + derivatives) is a degree 2n-1 method, so +it must reproduce t^3 from 3 nodes EXACTLY -- extrapolating to 0 gives 0, to roundoff. The old +Horner reconstruction walked the doubled node list at half speed, evaluating a polynomial that was +not the interpolant: this data came back ~1.6e-5 from 0, five orders above the samples' support. +*/ +BOOST_AUTO_TEST_CASE( hermite_reproduces_low_degree_polynomials_exactly ) +{ + DefaultPrecision(ambient_precision); + + BCT target_time(0,0); + unsigned int num_samples = 3; + + bertini::TimeCont times; + bertini::SampCont samples, derivatives; + Vec sample(1), derivative(1); + + for (auto const& t_str : {".1", ".05", ".025"}) + { + BCT t = ComplexFromString(t_str); + times.push_back(t); + sample << pow(t,3); samples.push_back(sample); // f(t) = t^3 + derivative << BCT(3)*pow(t,2); derivatives.push_back(derivative); // f'(t) = 3t^2 + } + + Vec approx = HermiteInterpolateAndSolve(target_time, num_samples, times, samples, derivatives); + BOOST_CHECK( norm(approx(0)) < 1e-20 ); +}//end hermite_reproduces_low_degree_polynomials_exactly + + @@ -217,9 +250,13 @@ BOOST_AUTO_TEST_CASE(hermite_interpolation) Vec< BCT > third_approx = HermiteInterpolateAndSolve(target_time,num_samples,times,samples,derivatives); - BOOST_CHECK((first_approx - correct).norm() < 1e-10); - BOOST_CHECK((second_approx - correct).norm() < 1e-10); - BOOST_CHECK((third_approx - correct).norm() < 1e-10); + // Tolerances calibrated to the TRUE Hermite interpolation errors of these windows (exact + // rational arithmetic): 1.162e-9, 4.539e-12, 1.773e-14 -- shrinking ~256x per halving. + // The old flat 1e-10 was calibrated to a mis-indexed Horner that happened to land closer + // to the truth than the actual interpolant does on the first window. + BOOST_CHECK((first_approx - correct).norm() < 2e-9); + BOOST_CHECK((second_approx - correct).norm() < 1e-11); + BOOST_CHECK((third_approx - correct).norm() < 1e-13); }//end hermite test case