fix(pseg): cycle-number bound is a ceiling; Hermite interpolation evaluates the actual interpolant

core/include/bertini2/endgames/interpolation.hpp

@@ -112,20 +112,19 @@ template<typename ComplexT>
 		}
 	}
 
-	//Start of Result from Hermite polynomial, this is using the diagonal of the 
-	//finite difference matrix.
-	Vec<ComplexT> 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<k}(x - z_j),  with a_k the
+	//divided-difference diagonal and z the DOUBLED node list (each time appears twice).
+	//Horner from the top: multiply by (x - z_k) before adding a_k, for k = 2n-2 down to 0.
+	//
+	//This loop previously indexed the doubled node list at half speed (z_ii paired with
+	//a_{2*ii}), evaluating a polynomial that is NOT the Hermite interpolant: it failed to
+	//reproduce even a cubic exactly.  The limit as the sample window slides to the target was
+	//unaffected, so the endgame still converged -- but at a degraded order, costing extra
+	//iterations and delivering poorer approximations than the samples support.
+	Vec<ComplexT> Result = space_differences(2*num_sample_points - 1,2*num_sample_points - 1);
+	for (int k = 2*static_cast<int>(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

core/include/bertini2/endgames/powerseries.hpp

@@ -378,13 +378,22 @@ class PowerSeriesEndgame :
 
 		RealT estimate = log(static_cast<RealT>(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<PowerSeriesConfig>();
+		if (estimate < 1)
 		  	upper_bound_on_cycle_number_ = 1;
 		else
 		{
-			using std::max;
-			auto upper_bound = unsigned(round(estimate)*this->template Get<PowerSeriesConfig>().cycle_number_amplification);
-			upper_bound_on_cycle_number_ = max(upper_bound,this->template Get<PowerSeriesConfig>().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<RealT>(ps_config.cycle_number_amplification);
+			if (!(amplified < static_cast<RealT>(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<unsigned>(amplified));
 		}
 
 		return upper_bound_on_cycle_number_;

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

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<BCT> times;
+	bertini::SampCont<BCT> samples, derivatives;
+	Vec<BCT> 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<BCT> 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
 
@@ -298,7 +335,7 @@ BOOST_AUTO_TEST_CASE(compute_bound_on_cycle_num)
 	my_endgame.ComputeBoundOnCycleNumber<BCT>();
 
 
-	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 +356,75 @@ BOOST_AUTO_TEST_CASE(compute_bound_on_cycle_num)
 	my_endgame.ComputeBoundOnCycleNumber<BCT>();
 
 
-	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
+
+
+/**
+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);
 
-} // end compute bound on cycle number 
+	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<BCT> times;
+	bertini::SampCont<BCT> samples;
+	Vec<BCT> 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<BCT>(1);
+
+	my_endgame.ComputeBoundOnCycleNumber<BCT>();
+	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<BCT>();
+	BOOST_CHECK_EQUAL(my_endgame.UpperBoundOnCycleNumber(), ceiling);
+} // end cycle_number_upper_bound_capped_for_near_unity_sample_ratios