diff --git a/test/test_1F1.hpp b/test/test_1F1.hpp index 049c9b06f..50903ed1e 100644 --- a/test/test_1F1.hpp +++ b/test/test_1F1.hpp @@ -160,7 +160,7 @@ void test_spots5(T, const char* type_name) template void test_spots6(T, const char* type_name) { - static const boost::array, 75> hypergeometric_1F1_bugs = { { + static const boost::array, 91> hypergeometric_1F1_bugs = { { { { static_cast(17955.561660766602), static_cast(9.6968994205831605e-09), static_cast(-82.406154185533524), SC_(6.98056008378736714088730927132364938220428678e-11) }}, { { static_cast(17955.561660766602), static_cast(-9.6968994205831605e-09), static_cast(-82.406154185533524), SC_(-6.98055306629610746072607353939306734740549551e-11) }}, { { static_cast(-17955.561660766602), static_cast(-9.6968994205831605e-09), static_cast(82.406154185533524), SC_(-42897094853118832762870100.8669248353530950866) }} , @@ -267,6 +267,24 @@ void test_spots6(T, const char* type_name) { { std::ldexp((double)10485122560373760, -40), std::ldexp((double)-11098279821997376, -39), std::ldexp((double)16925852270592000, -45), SC_(9.77378642649349178995585980824930703376759021e-98) } }, { { std::ldexp((double)10485292967829248, -40), std::ldexp((double)-14859721380002656, -35), std::ldexp((double)13729956970496000, -44), SC_(3.41094899910311302761937103011397882987669395e-08) } }, { { std::ldexp((double)10485037389193216, -40), std::ldexp((double)-10840488483391544, -35), std::ldexp((double)17577061875712000, -45), SC_(2.8030884395368690164859926372380406504460219e-07) } }, + // + // Negative a and b worst cases: + { { std::ldexp((double)-9281686323200000, -44), std::ldexp((double)-14062138056704000, -44), std::ldexp((double)13563284652032000, -44), SC_(2.83381029611748905436564183928049453625733389e+265) } }, + { { std::ldexp((double)-17049048150016000, -44), std::ldexp((double)-16971363917824000, -45), std::ldexp((double)11759598960640000, -49), SC_(4636596575297708282.15391199523219256661832333) }}, + { { std::ldexp((double)-14233964060672000, -45), std::ldexp((double)-12648356216832000, -47), std::ldexp((double)9597206757376000, -46), SC_(-1.29952965544474451915331906705211320124261355e+104) }}, + { { std::ldexp((double)-16705334214656000, -45), std::ldexp((double)-15447756718080000, -46), std::ldexp((double)16395884134400000, -47), SC_(5.40680141346616359293013198457674276017117348e+113) }}, + { { std::ldexp((double)-13991530405888000, -45), std::ldexp((double)-10196587347968000, -46), std::ldexp((double)13331347734528000, -46), SC_(1.28612752976615347819085089716937824474111365e+138) }}, + { { std::ldexp((double)-15134950760448000, -45), std::ldexp((double)-14587193786368000, -48), std::ldexp((double)17022855921664000, -46), SC_(-8.81689040877580073465185463207591010592963947e+115) }}, + { { std::ldexp((double)-14854672039936000, -45), std::ldexp((double)-10436558200832000, -45), std::ldexp((double)11370918969344000, -47), SC_(8.85535247272534117445528460568920680343270651e+54) } }, + { { std::ldexp((double)-16711069286400000, -46), std::ldexp((double)-14809815056384000, -46), std::ldexp((double)10469312954368000, -47), SC_(50343352353398198766339890377962553344.0) } }, + { { std::ldexp((double)-15026786402304000, -45), std::ldexp((double)-16687356968960000, -46), std::ldexp((double)14895621603328000, -47), SC_(2.85329560424602656900599696665580727040444839e+95) } }, + { { std::ldexp((double)-15519073435648000, -45), std::ldexp((double)-14162009718784000, -45), std::ldexp((double)9997818855424000, -48), SC_(95767987018108517.7639995774279100260173436254) } }, + { { std::ldexp((double)-15317481275392000, -46), std::ldexp((double)-16531865931776000, -44), std::ldexp((double)17586268880896000, -45), SC_(1.47012480470837240834700634331053708005164679e+104) }}, + { { std::ldexp((double)-11335669673984000, -44), std::ldexp((double)-13146047094784000, -44), std::ldexp((double)13671437864960000, -44), SC_(-2.18876072849870893180776779453843962111085257e+288) }}, + { { std::ldexp((double)-16877985234944000, -46), std::ldexp((double)-14384006086656000, -46), std::ldexp((double)9074349342720000, -47), SC_(15376193613462463541358751744655360.0) }}, + { { std::ldexp((double)-9751199809536000, -45), std::ldexp((double)-17654191685632000, -47), std::ldexp((double)10587451850752000, -47), SC_(-1.9601415510439595625538337964298353914980331e+68) }}, + { { std::ldexp((double)-15233620754432000, -45), std::ldexp((double)-12708283072512000, -46), std::ldexp((double)10255461007360000, -46), SC_(-5.43441063616790758618595678580168066778261756e+125) }}, + { { std::ldexp((double)-11241354149888000, -45), std::ldexp((double)-9580579905536000, -45), std::ldexp((double)12224976846848000, -47), SC_(1.20468565484700674058707264904280026536015299e+46) }}, } }; static const boost::array, 2> hypergeometric_1F1_big_bugs = { { #if DBL_MAX_EXP == LDBL_MAX_EXP diff --git a/tools/hypergeometric_1F1_error_plot.cpp b/tools/hypergeometric_1F1_error_plot.cpp index 497a003b1..649cd7bed 100644 --- a/tools/hypergeometric_1F1_error_plot.cpp +++ b/tools/hypergeometric_1F1_error_plot.cpp @@ -122,7 +122,7 @@ int main() mpfr_float mp_expected; try { - mp_expected = boost::math::hypergeometric_pFq_precision({ mpfr_float(a) }, { mpfr_float(b) }, mpfr_float(z), 25, 3.0); + mp_expected = boost::math::hypergeometric_pFq_precision({ mpfr_float(a) }, { mpfr_float(b) }, mpfr_float(z), 25, 200.0); expected = (test_type)mp_expected; } catch (const std::exception&) diff --git a/tools/hypergeometric_1F1_map_neg_b_fwd_recurrence.cpp b/tools/hypergeometric_1F1_map_neg_b_fwd_recurrence.cpp index cfccf0dea..9655ddd57 100644 --- a/tools/hypergeometric_1F1_map_neg_b_fwd_recurrence.cpp +++ b/tools/hypergeometric_1F1_map_neg_b_fwd_recurrence.cpp @@ -56,7 +56,7 @@ int main() std::cin >> b_mult; double error_limit = 200; - double time_limit = 100.0; + double time_limit = 10.0; for (test_type a = a_start; a < a_end; a_start < 0 ? a /= a_mult : a *= a_mult) { @@ -65,9 +65,10 @@ int main() test_type z_mult = 2; test_type last_good = 0; test_type bad = 0; - for (test_type z = 1; z < 1e10; z *= z_mult, z_mult *= 2) - { - try { + try { + for (test_type z = 1; z < 1e10; z *= z_mult, z_mult *= 2) + { + // std::cout << "z = " << z << std::endl; boost::uintmax_t max_iter = 1000; test_type calc = boost::math::tools::function_ratio_from_forwards_recurrence(boost::math::detail::hypergeometric_1F1_recurrence_a_and_b_coefficients(a, b, z), std::numeric_limits::epsilon() * 2, max_iter); test_type reference = (test_type)(boost::math::hypergeometric_pFq_precision({ mpfr_float(a) }, { mpfr_float(b) }, mpfr_float(z), 50, time_limit) / boost::math::hypergeometric_pFq_precision({ mpfr_float(a + 1) }, { mpfr_float(b + 1) }, mpfr_float(z), std::numeric_limits::digits10 * 2, time_limit)); @@ -83,11 +84,11 @@ int main() bad = z; } } - catch (const std::exception& e) - { - std::cout << "Unexpected exception: " << e.what() << std::endl; - std::cout << "For a = " << a << " b = " << b << " z = " << z << std::endl; - } + } + catch (const std::exception& e) + { + std::cout << "Unexpected exception: " << e.what() << std::endl; + std::cout << "For a = " << a << " b = " << b << " z = " << bad * z_mult / 2 << std::endl; } test_type z_limit; if (0 == bad) @@ -109,14 +110,16 @@ int main() }, bad, last_good, boost::math::tools::equal_floor()).first; z_limit = floor(z_limit + 2); // Give ourselves some headroom! } + // std::cout << "z_limit = " << z_limit << std::endl; // // Now over again for backwards recurrence domain at the same points: // - bad = z_limit; + bad = z_limit > 1e10 ? 1e10 : z_limit; last_good = 0; z_mult = 1.1; - for (test_type z = z_limit; z > 1; z /= z_mult, z_mult *= 2) + for (test_type z = bad; z > 1; z /= z_mult, z_mult *= 2) { + // std::cout << "z = " << z << std::endl; try { boost::uintmax_t max_iter = 1000; test_type calc = boost::math::tools::function_ratio_from_backwards_recurrence(boost::math::detail::hypergeometric_1F1_recurrence_a_and_b_coefficients(a, b, z), std::numeric_limits::epsilon() * 2, max_iter); @@ -135,6 +138,7 @@ int main() } catch (const std::exception& e) { + bad = z; std::cout << "Unexpected exception: " << e.what() << std::endl; std::cout << "For a = " << a << " b = " << b << " z = " << z << std::endl; }