In perl.git, the branch blead has been updated <http://perl5.git.perl.org/perl.git/commitdiff/e64e4e046f94ead9bf3ab016f056daf1e01ec312?hp=51a0624d13d30bf656b67b2f094c54e7ab0fd859>
- Log ----------------------------------------------------------------- commit e64e4e046f94ead9bf3ab016f056daf1e01ec312 Author: Jarkko Hietaniemi <[email protected]> Date: Tue Sep 2 20:44:21 2014 -0400 POSIX math: add tests. M ext/POSIX/t/math.t commit d867902f8a7f25393479a46177c69edb59747096 Author: Jarkko Hietaniemi <[email protected]> Date: Tue Sep 2 20:33:53 2014 -0400 POSIX math: add isinf/isnan/isfinite tests. M ext/POSIX/t/math.t commit ffbadafbbdde319912a03812571f9f5a17e57853 Author: Jarkko Hietaniemi <[email protected]> Date: Tue Sep 2 19:47:11 2014 -0400 POSIX math: erfc implementation. (Chiani, Dardari & Simon (2003), via Wikipedia) M ext/POSIX/POSIX.xs M ext/POSIX/t/math.t commit cc06cdb1aa1976f96fe7892fa1691277d8a16543 Author: Jarkko Hietaniemi <[email protected]> Date: Tue Sep 2 19:36:35 2014 -0400 POSIX math: looser error bound for erf test. M ext/POSIX/t/math.t commit 6727dab6a8ce835e67a5670f1ebbce0c43e589bf Author: Jarkko Hietaniemi <[email protected]> Date: Tue Sep 2 19:26:36 2014 -0400 POSIX math: more hypot tests. M ext/POSIX/t/math.t commit 3e2e323f162f541949f6daec0adaca39172b0a8d Author: Jarkko Hietaniemi <[email protected]> Date: Tue Sep 2 19:21:24 2014 -0400 POSIX math: isless etc tests. M ext/POSIX/t/math.t commit 7bbf2c9b18e21dae19461e18dc4f9f63211a6fc3 Author: Jarkko Hietaniemi <[email protected]> Date: Tue Sep 2 19:17:41 2014 -0400 POSIX math: better hypot. M ext/POSIX/POSIX.xs commit 7895eac5c8bcf075dcdc364b3b43be094c809502 Author: Jarkko Hietaniemi <[email protected]> Date: Tue Sep 2 19:10:36 2014 -0400 POSIX math: erf implementation from http://www.johndcook.com/ (Also add sources to other formulas where they weren't obvious.) M ext/POSIX/POSIX.xs M ext/POSIX/t/math.t commit 26e0f0d373d8649f6ddc462743c555bd39166946 Author: Jarkko Hietaniemi <[email protected]> Date: Tue Sep 2 17:48:25 2014 -0400 POSIX math: Define my_fpclassify only if no c99_classify already. ("c99_classify already" = we saw HAS_FPCLASSIFY and fpclassify macro, which together smell like C99) M ext/POSIX/POSIX.xs ----------------------------------------------------------------------- Summary of changes: ext/POSIX/POSIX.xs | 72 +++++++++++++++++++++++++++++++++++++++++++++--------- ext/POSIX/t/math.t | 33 ++++++++++++++++++++----- 2 files changed, 87 insertions(+), 18 deletions(-) diff --git a/ext/POSIX/POSIX.xs b/ext/POSIX/POSIX.xs index 77d4ec8..2b42a17 100644 --- a/ext/POSIX/POSIX.xs +++ b/ext/POSIX/POSIX.xs @@ -139,17 +139,17 @@ log1p log2 logb lrint nan nearbyint nextafter nexttoward remainder remquo rint round scalbn signbit tgamma trunc - * Berkeley/SVID extensions: + * Berkeley/SVID extensions: - j0 j1 jn y0 y1 yn + j0 j1 jn y0 y1 yn - * Configure already (5.21.0) scans for: + * Configure already (5.21.0) scans for: - fpclassify isfinite isinf isnan ilogb*l* signbit + fpclassify isfinite isinf isnan ilogb*l* signbit - * For floating-point round mode (which matters for e.g. lrint and rint) + * For floating-point round mode (which matters for e.g. lrint and rint) - fegetround fesetround + fegetround fesetround */ @@ -464,8 +464,38 @@ static NV my_copysign(NV x, NV y) /* XXX cosh (though c89) */ -/* XXX erf -- non-trivial */ -/* XXX erfc -- non-trivial */ +#ifndef c99_erf +static NV my_erf(NV x) +{ + /* http://www.johndcook.com/cpp_erf.html -- public domain */ + NV a1 = 0.254829592; + NV a2 = -0.284496736; + NV a3 = 1.421413741; + NV a4 = -1.453152027; + NV a5 = 1.061405429; + NV p = 0.3275911; + + int sign = x < 0 ? -1 : 1; /* Save the sign. */ + x = PERL_ABS(x); + + /* Abramowitz and Stegun formula 7.1.26 */ + NV t = 1.0 / (1.0 + p * x); + NV y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1) * t * exp(-x*x); + + return sign * y; +} +# define c99_erf my_erf +#endif + +/* While in theory erfc(x) is just 1-erf(x), thanks to numerical stability + * things are not so easy. */ +#ifndef c99_erfc +static NV my_erfc(NV x) { + /* Chiani, Dardari & Simon (2003), via Wikipedia */ + return Perl_exp(-x*x) / 6.0 + Perl_exp(-4.0/3.0 * x*x) / 2.0; +} +# define c99_erfc my_erfc +#endif #ifndef c99_exp2 static NV my_exp2(NV x) @@ -479,6 +509,8 @@ static NV my_exp2(NV x) static NV my_expm1(NV x) { if (PERL_ABS(x) < 1e-5) + /* http://www.johndcook.com/cpp_expm1.html -- public domain. + * Also including the cubic term. */ /* Probably not enough for long doubles. */ return x * (1.0 + x * (0.5 + x / 6.0)); /* Taylor series */ else @@ -521,6 +553,8 @@ static NV my_fmin(NV x, NV y) # define c99_fmin my_fmin #endif +#ifndef c99_fpclassify + static IV my_fpclassify(NV x) { #if defined(HAS_FPCLASSIFY) && defined(FP_PLUS_INF) /* E.g. HP-UX */ @@ -590,14 +624,26 @@ static IV my_fpclassify(NV x) #endif } +#endif + #ifndef c99_hypot static NV my_hypot(NV x, NV y) { - if (x > 0.0) { - NV t = y / x; - return PERL_ABS(x) * Perl_sqrt(1 + t * t); + /* http://en.wikipedia.org/wiki/Hypot */ + NV t; + x = PERL_ABS(x); /* Take absolute values. */ + if (y == 0) + return x; + if (Perl_isnan(y)) + return NV_INF; + y = PERL_ABS(y); + if (x < y) { /* Swap so that y is less. */ + t = x; + x = y; + y = t; } - return NV_NAN; + t = y / x; + return x * sqrt(1.0 + t * t); } # define c99_hypot my_hypot #endif @@ -615,6 +661,8 @@ static IV my_ilogb(NV x) #ifndef c99_log1p static NV my_log1p(NV x) { + /* http://www.johndcook.com/cpp_log_one_plus_x.html -- public domain. + * Including also quadratic term. */ if (PERL_ABS(x) > 1e-4) return Perl_log(1.0 + x); else diff --git a/ext/POSIX/t/math.t b/ext/POSIX/t/math.t index 0f0642b..043a578 100644 --- a/ext/POSIX/t/math.t +++ b/ext/POSIX/t/math.t @@ -64,7 +64,9 @@ SKIP: { cmp_ok(abs(M_SQRT2 - 1.4142135623731), '<', 1e-9, "M_SQRT2"); cmp_ok(abs(M_E - 2.71828182845905), '<', 1e-9, "M_E"); cmp_ok(abs(M_PI - 3.14159265358979), '<', 1e-9, "M_PI"); + cmp_ok(abs(acosh(2) - 1.31695789692482), '<', 1e-9, "acosh"); cmp_ok(abs(asinh(1) - 0.881373587019543), '<', 1e-9, "asinh"); + cmp_ok(abs(atanh(0.5) - 0.549306144334055), '<', 1e-9, "atanh"); cmp_ok(abs(cbrt(8) - 2), '<', 1e-9, "cbrt"); cmp_ok(abs(cbrt(-27) - -3), '<', 1e-9, "cbrt"); is(copysign(3.14, -2), -3.14, "copysign"); @@ -84,26 +86,36 @@ SKIP: { is(fpclassify(NAN), FP_NAN, "fpclassify NAN"); } is(hypot(3, 4), 5, "hypot 3 4"); + is(hypot(-2, 1), sqrt(5), "hypot -1 2"); is(ilogb(255), 7, "ilogb 255"); is(ilogb(256), 8, "ilogb 256"); SKIP: { unless ($Config{d_isfinite}) { - skip "no isfinite", 1; + skip "no isfinite", 3; } - ok(isfinite(1), "isfinite"); + ok(isfinite(1), "isfinite 1"); + ok(!isfinite(Inf), "isfinite Inf"); + ok(!isfinite(NaN), "isfinite NaN"); } SKIP: { unless ($Config{d_isinf}) { - skip "no isinf", 1; + skip "no isinf", 4; } - ok(isinf(INFINITY), "isinf"); + ok(isinf(INFINITY), "isinf INFINITY"); + ok(isinf(Inf), "isinf Inf"); + ok(!isinf(NaN), "isinf NaN"); + ok(!isinf(42), "isinf 42"); } SKIP: { unless ($Config{d_isnan}) { - skip "no isnan", 1; + skip "no isnan", 4; } - ok(isnan(NAN), "isnan"); + ok(isnan(NAN), "isnan NAN"); + ok(isnan(NaN), "isnan NaN"); + ok(!isnan(Inf), "isnan Inf"); + ok(!isnan(42), "isnan Inf"); } + cmp_ok(nan(), '!=', nan(), 'nan'); cmp_ok(abs(log1p(2) - 1.09861228866811), '<', 1e-9, "log1p"); cmp_ok(abs(log1p(1e-6) - 9.99999500000333e-07), '<', 1e-9, "log1p"); cmp_ok(abs(log2(8) - 3), '<', 1e-9, "log2"); @@ -126,6 +138,15 @@ SKIP: { is(trunc(-2.5), -2, "trunc -2.5"); is(trunc(2.75), 2, "trunc 2.75"); is(trunc(-2.75), -2, "trunc -2.75"); + ok(isless(1, 2), "isless 1 2"); + ok(!isless(2, 1), "isless 2 1"); + ok(!isless(1, 1), "isless 1 1"); + ok(!isless(1, NaN), "isless 1 NaN"); + ok(isgreater(2, 1), "isgreater 2 1"); + ok(islessequal(1, 1), "islessequal 1 1"); + ok(isunordered(1, NaN), "isunordered 1 NaN"); + cmp_ok(abs(erf(1) - 0.842700792949715), '<', 1.5e-7, "erf 1"); + cmp_ok(abs(erfc(1) - 0.157299207050285), '<', 1.5e-7, "erfc 1"); } done_testing(); -- Perl5 Master Repository
