diff --git a/src/distribution/log_normal.rs b/src/distribution/log_normal.rs index d36dc504..7628c589 100644 --- a/src/distribution/log_normal.rs +++ b/src/distribution/log_normal.rs @@ -750,13 +750,13 @@ mod tests { test_absolute(-0.1, 0.1, 1.0, 1e-107, sf(0.1)); // Wolfram Alpha:: SurvivalFunction[ LogNormalDistribution(-0.1, 0.1), 0.8] - test_absolute(-0.1, 0.1, 0.890919989231123, 1e-14, sf(0.8)); + test_absolute(-0.1, 0.1, 0.890919989236242017902656326193, 1e-14, sf(0.8)); // Wolfram Alpha:: SurvivalFunction[LogNormalDistribution[1.5, 1], 0.8] - test_absolute(1.5, 1.0, 0.957568715612642, 1e-14, sf(0.8)); + test_absolute(1.5, 1.0, 0.957568715614411289470712253063, 1e-14, sf(0.8)); // Wolfram Alpha:: SurvivalFunction[ LogNormalDistribution(2.5, 1.5), 0.1] - test_absolute(2.5, 1.5, 0.9993169594777358, 1e-14, sf(0.1)); + test_absolute(2.5, 1.5, 0.999316959477792115067727898746, 1e-14, sf(0.1)); } #[test] diff --git a/src/function/erf.rs b/src/function/erf.rs index 63238712..023670df 100644 --- a/src/function/erf.rs +++ b/src/function/erf.rs @@ -594,82 +594,94 @@ fn erf_impl(z: f64, inv: bool) -> f64 { ( evaluate::polynomial(z - 0.5, ERF_IMPL_BN) / evaluate::polynomial(z - 0.5, ERF_IMPL_BD), - 0.3440242112, + 0.3440242111682891845703125, ) } else if z < 1.25 { ( evaluate::polynomial(z - 0.75, ERF_IMPL_CN) / evaluate::polynomial(z - 0.75, ERF_IMPL_CD), - 0.419990927, + 0.4199909269809722900390625, ) } else if z < 2.25 { ( evaluate::polynomial(z - 1.25, ERF_IMPL_DN) / evaluate::polynomial(z - 1.25, ERF_IMPL_DD), - 0.4898625016, + 0.489862501621246337890625, ) } else if z < 3.5 { ( evaluate::polynomial(z - 2.25, ERF_IMPL_EN) / evaluate::polynomial(z - 2.25, ERF_IMPL_ED), - 0.5317370892, + 0.5317370891571044921875, ) } else if z < 5.25 { ( evaluate::polynomial(z - 3.5, ERF_IMPL_FN) / evaluate::polynomial(z - 3.5, ERF_IMPL_FD), - 0.5489973426, + 0.548997342586517333984375, ) } else if z < 8.0 { ( evaluate::polynomial(z - 5.25, ERF_IMPL_GN) / evaluate::polynomial(z - 5.25, ERF_IMPL_GD), - 0.5571740866, + 0.55717408657073974609375, ) } else if z < 11.5 { ( evaluate::polynomial(z - 8.0, ERF_IMPL_HN) / evaluate::polynomial(z - 8.0, ERF_IMPL_HD), - 0.5609807968, + 0.56098079681396484375, ) } else if z < 17.0 { ( evaluate::polynomial(z - 11.5, ERF_IMPL_IN) / evaluate::polynomial(z - 11.5, ERF_IMPL_ID), - 0.5626493692, + 0.56264936923980712890625, ) } else if z < 24.0 { ( evaluate::polynomial(z - 17.0, ERF_IMPL_JN) / evaluate::polynomial(z - 17.0, ERF_IMPL_JD), - 0.5634598136, + 0.563459813594818115234375, ) } else if z < 38.0 { ( evaluate::polynomial(z - 24.0, ERF_IMPL_KN) / evaluate::polynomial(z - 24.0, ERF_IMPL_KD), - 0.5638477802, + 0.5638477802276611328125, ) } else if z < 60.0 { ( evaluate::polynomial(z - 38.0, ERF_IMPL_LN) / evaluate::polynomial(z - 38.0, ERF_IMPL_LD), - 0.5640528202, + 0.5640528202056884765625, ) } else if z < 85.0 { ( evaluate::polynomial(z - 60.0, ERF_IMPL_MN) / evaluate::polynomial(z - 60.0, ERF_IMPL_MD), - 0.5641309023, + 0.56413090229034423828125, ) } else { ( evaluate::polynomial(z - 85.0, ERF_IMPL_NN) / evaluate::polynomial(z - 85.0, ERF_IMPL_ND), - 0.5641584396, + 0.56415843963623046875, ) }; - let g = (-z * z).exp() / z; + // `z * z` is rounded before `exp` sees it, which costs roughly `z^2` + // ulps in the result (~460 ulps by z = 27). Recover the rounding error + // of the square exactly with a Dekker product and fold it back in: + // `exp(-z^2) = exp(-sq) * exp(-err)`, and `exp(-err) ~= 1 - err` since + // `|err| <= ulp(sq)/2`. Dekker's split is used instead of + // `f64::mul_add` because the latter falls back to a slow software FMA + // on targets without the hardware instruction (e.g. baseline x86-64). + let sq = z * z; + let split = 134_217_729.0 * z; // 2^27 + 1; cannot overflow, z < 110 + let z_hi = split - (split - z); + let z_lo = z - z_hi; + let err = ((z_hi * z_hi - sq) + 2.0 * z_hi * z_lo) + z_lo * z_lo; + let g = (-sq).exp() * (1.0 - err) / z; g * b + g * r } else { 0.0 @@ -747,19 +759,19 @@ mod tests { #[test] fn test_erf() { assert!(erf(f64::NAN).is_nan()); - prec::assert_abs_diff_eq!(erf(-1.0), -0.84270079294971486934122063508260925929606699796630291, epsilon = 1e-11); + prec::assert_abs_diff_eq!(erf(-1.0), -0.84270079294971486934122063508260925929606699796630291, epsilon = 4.441e-16); assert_eq!(erf(0.0), 0.0); assert_eq!(erf(1e-15), 0.0000000000000011283791670955126615773132947717431253912942469337536); assert_eq!(erf(0.1), 0.1124629160182848984047122510143040617233925185058162); prec::assert_abs_diff_eq!(erf(0.2), 0.22270258921047846617645303120925671669511570710081967, epsilon = 1e-16); assert_eq!(erf(0.3), 0.32862675945912741618961798531820303325847175931290341); assert_eq!(erf(0.4), 0.42839235504666847645410962730772853743532927705981257); - prec::assert_abs_diff_eq!(erf(0.5), 0.5204998778130465376827466538919645287364515757579637, epsilon = 1e-9); - prec::assert_abs_diff_eq!(erf(1.0), 0.84270079294971486934122063508260925929606699796630291, epsilon = 1e-11); - prec::assert_abs_diff_eq!(erf(1.5), 0.96610514647531072706697626164594785868141047925763678, epsilon = 1e-11); - prec::assert_abs_diff_eq!(erf(2.0), 0.99532226501895273416206925636725292861089179704006008, epsilon = 1e-11); - prec::assert_abs_diff_eq!(erf(2.5), 0.99959304798255504106043578426002508727965132259628658, epsilon = 1e-13); - prec::assert_abs_diff_eq!(erf(3.0), 0.99997790950300141455862722387041767962015229291260075, epsilon = 1e-11); + prec::assert_abs_diff_eq!(erf(0.5), 0.5204998778130465376827466538919645287364515757579637, epsilon = 4.441e-16); + prec::assert_abs_diff_eq!(erf(1.0), 0.84270079294971486934122063508260925929606699796630291, epsilon = 4.441e-16); + prec::assert_abs_diff_eq!(erf(1.5), 0.96610514647531072706697626164594785868141047925763678, epsilon = 4.441e-16); + prec::assert_abs_diff_eq!(erf(2.0), 0.99532226501895273416206925636725292861089179704006008, epsilon = 4.441e-16); + prec::assert_abs_diff_eq!(erf(2.5), 0.99959304798255504106043578426002508727965132259628658, epsilon = 4.441e-16); + prec::assert_abs_diff_eq!(erf(3.0), 0.99997790950300141455862722387041767962015229291260075, epsilon = 4.441e-16); assert_eq!(erf(4.0), 0.99999998458274209971998114784032651311595142785474641); assert_eq!(erf(5.0), 0.99999999999846254020557196514981165651461662110988195); assert_eq!(erf(6.0), 0.99999999999999997848026328750108688340664960081261537); @@ -770,24 +782,24 @@ mod tests { #[test] fn test_erfc() { assert!(erfc(f64::NAN).is_nan()); - prec::assert_abs_diff_eq!(erfc(-1.0), 1.8427007929497148693412206350826092592960669979663028, epsilon = 1e-11); + prec::assert_abs_diff_eq!(erfc(-1.0), 1.8427007929497148693412206350826092592960669979663028, epsilon = 8.882e-16); assert_eq!(erfc(0.0), 1.0); - prec::assert_abs_diff_eq!(erfc(0.1), 0.88753708398171510159528774898569593827660748149418343, epsilon = 1e-15); + prec::assert_abs_diff_eq!(erfc(0.1), 0.88753708398171510159528774898569593827660748149418343, epsilon = 4.441e-16); assert_eq!(erfc(0.2), 0.77729741078952153382354696879074328330488429289918085); assert_eq!(erfc(0.3), 0.67137324054087258381038201468179696674152824068709621); - prec::assert_abs_diff_eq!(erfc(0.4), 0.57160764495333152354589037269227146256467072294018715, epsilon = 1e-15); - prec::assert_abs_diff_eq!(erfc(0.5), 0.47950012218695346231725334610803547126354842424203654, epsilon = 1e-9); - prec::assert_abs_diff_eq!(erfc(1.0), 0.15729920705028513065877936491739074070393300203369719, epsilon = 1e-11); - prec::assert_abs_diff_eq!(erfc(1.5), 0.033894853524689272933023738354052141318589520742363247, epsilon = 1e-11); - prec::assert_abs_diff_eq!(erfc(2.0), 0.0046777349810472658379307436327470713891082029599399245, epsilon = 1e-11); - prec::assert_abs_diff_eq!(erfc(2.5), 0.00040695201744495893956421573997491272034867740371342016, epsilon = 1e-13); - prec::assert_abs_diff_eq!(erfc(3.0), 0.00002209049699858544137277612958232037984770708739924966, epsilon = 1e-11); - prec::assert_abs_diff_eq!(erfc(4.0), 0.000000015417257900280018852159673486884048572145253589191167, epsilon = 1e-18); - prec::assert_abs_diff_eq!(erfc(5.0), 0.0000000000015374597944280348501883434853833788901180503147233804, epsilon = 1e-22); - prec::assert_abs_diff_eq!(erfc(6.0), 2.1519736712498913116593350399187384630477514061688559e-17, epsilon = 1e-26); - prec::assert_abs_diff_eq!(erfc(10.0), 2.0884875837625447570007862949577886115608181193211634e-45, epsilon = 1e-55); - prec::assert_abs_diff_eq!(erfc(15.0), 7.2129941724512066665650665586929271099340909298253858e-100, epsilon = 1e-109); - prec::assert_abs_diff_eq!(erfc(20.0), 5.3958656116079009289349991679053456040882726709236071e-176, epsilon = 1e-186); + prec::assert_abs_diff_eq!(erfc(0.4), 0.57160764495333152354589037269227146256467072294018715, epsilon = 4.441e-16); + prec::assert_abs_diff_eq!(erfc(0.5), 0.47950012218695346231725334610803547126354842424203654, epsilon = 2.220e-16); + prec::assert_abs_diff_eq!(erfc(1.0), 0.15729920705028513065877936491739074070393300203369719, epsilon = 1.110e-16); + prec::assert_abs_diff_eq!(erfc(1.5), 0.033894853524689272933023738354052141318589520742363247, epsilon = 2.776e-17); + prec::assert_abs_diff_eq!(erfc(2.0), 0.0046777349810472658379307436327470713891082029599399245, epsilon = 3.469e-18); + prec::assert_abs_diff_eq!(erfc(2.5), 0.00040695201744495893956421573997491272034867740371342016, epsilon = 2.168e-19); + prec::assert_abs_diff_eq!(erfc(3.0), 0.00002209049699858544137277612958232037984770708739924966, epsilon = 1.355e-20); + prec::assert_abs_diff_eq!(erfc(4.0), 0.000000015417257900280018852159673486884048572145253589191167, epsilon = 1.323e-23); + prec::assert_abs_diff_eq!(erfc(5.0), 0.0000000000015374597944280348501883434853833788901180503147233804, epsilon = 8.078e-28); + prec::assert_abs_diff_eq!(erfc(6.0), 2.1519736712498913116593350399187384630477514061688559e-17, epsilon = 1.233e-32); + prec::assert_abs_diff_eq!(erfc(10.0), 2.0884875837625447570007862949577886115608181193211634e-45, epsilon = 1.245e-60); + prec::assert_abs_diff_eq!(erfc(15.0), 7.2129941724512066665650665586929271099340909298253858e-100, epsilon = 4.061e-115); + prec::assert_abs_diff_eq!(erfc(20.0), 5.3958656116079009289349991679053456040882726709236071e-176, epsilon = 2.806e-191); assert_eq!(erfc(30.0), 2.5646562037561116000333972775014471465488897227786155e-393); assert_eq!(erfc(50.0), 2.0709207788416560484484478751657887929322509209953988e-1088); assert_eq!(erfc(80.0), 2.3100265595063985852034904366341042118385080919280966e-2782); @@ -800,9 +812,9 @@ mod tests { assert!(erf_inv(f64::NAN).is_nan()); assert_eq!(erf_inv(-1.0), f64::NEG_INFINITY); assert_eq!(erf_inv(0.0), 0.0); - prec::assert_abs_diff_eq!(erf_inv(1e-15), 8.86226925452758013649e-16, epsilon = 1e-30); + prec::assert_abs_diff_eq!(erf_inv(1e-15), 8.86226925452758013649e-16, epsilon = 3.944e-31); assert_eq!(erf_inv(0.1), 0.08885599049425768701574); - prec::assert_abs_diff_eq!(erf_inv(0.2), 0.1791434546212916764927, epsilon = 1e-15); + prec::assert_abs_diff_eq!(erf_inv(0.2), 0.1791434546212916764927, epsilon = 1.110e-16); assert_eq!(erf_inv(0.3), 0.272462714726754355622); assert_eq!(erf_inv(0.4), 0.3708071585935579290582); assert_eq!(erf_inv(0.5), 0.4769362762044698733814); @@ -814,13 +826,13 @@ mod tests { #[test] fn test_erfc_inv() { assert_eq!(erfc_inv(0.0), f64::INFINITY); - prec::assert_abs_diff_eq!(erfc_inv(1e-100), 15.065574702593, epsilon = 1e-11); - prec::assert_abs_diff_eq!(erfc_inv(1e-30), 8.1486162231699, epsilon = 1e-12); - prec::assert_abs_diff_eq!(erfc_inv(1e-20), 6.6015806223551, epsilon = 1e-13); - prec::assert_abs_diff_eq!(erfc_inv(1e-10), 4.5728249585449249378479309946884581365517663258840893, epsilon = 1e-7); - prec::assert_abs_diff_eq!(erfc_inv(1e-5), 3.1234132743415708640270717579666062107939039971365252, epsilon = 1e-11); - prec::assert_abs_diff_eq!(erfc_inv(0.1), 1.1630871536766741628440954340547000483801487126688552, epsilon = 1e-14); - prec::assert_abs_diff_eq!(erfc_inv(0.2), 0.90619380243682330953597079527631536107443494091638384, epsilon = 1e-15); + prec::assert_abs_diff_eq!(erfc_inv(1e-100), 15.065574702592645704404610541369, epsilon = 7.105e-15); + prec::assert_abs_diff_eq!(erfc_inv(1e-30), 8.1486162231698646073845666606481, epsilon = 7.105e-15); + prec::assert_abs_diff_eq!(erfc_inv(1e-20), 6.6015806223551425615163916324187, epsilon = 3.553e-15); + prec::assert_abs_diff_eq!(erfc_inv(1e-10), 4.5728249673894852787410436731442, epsilon = 3.553e-15); + prec::assert_abs_diff_eq!(erfc_inv(1e-5), 3.1234132743408750302472925681804, epsilon = 1.776e-15); + prec::assert_abs_diff_eq!(erfc_inv(0.1), 1.1630871536766740867262542605629, epsilon = 8.882e-16); + prec::assert_abs_diff_eq!(erfc_inv(0.2), 0.90619380243682322007116270309566, epsilon = 4.441e-16); assert_eq!(erfc_inv(0.5), 0.47693627620446987338141835364313055980896974905947083); assert_eq!(erfc_inv(1.0), 0.0); assert_eq!(erfc_inv(1.5), -0.47693627620446987338141835364313055980896974905947083);