diff options
author | Jarkko Hietaniemi <jhi@iki.fi> | 2014-08-29 22:11:58 -0400 |
---|---|---|
committer | Jarkko Hietaniemi <jhi@iki.fi> | 2014-08-31 17:53:06 -0400 |
commit | 39b5f1c42e8171dd236a435af20ce1ff13c1af47 (patch) | |
tree | 949cf5d7c7f08b2d8cdcb22ec84d60947c5f8bb0 | |
parent | d2bce0fdf1a48c474f385b70ce22910b7b7a46a7 (diff) | |
download | perl-39b5f1c42e8171dd236a435af20ce1ff13c1af47.tar.gz |
POSIX math: Portability emulations and constants.
Plus fix the cmp_ok tests which had epsilon 18 magnitudes off.
(Didn't cause any false positives, luckily.)
-rw-r--r-- | ext/POSIX/POSIX.xs | 241 | ||||
-rw-r--r-- | ext/POSIX/t/math.t | 22 |
2 files changed, 251 insertions, 12 deletions
diff --git a/ext/POSIX/POSIX.xs b/ext/POSIX/POSIX.xs index fa9c8c6de0..aac26082e5 100644 --- a/ext/POSIX/POSIX.xs +++ b/ext/POSIX/POSIX.xs @@ -54,6 +54,46 @@ #include <unistd.h> #endif +#ifndef M_E +# define M_E 2.71828182845904523536028747135266250 +#endif +#ifndef M_LOG2E +# define M_LOG2E 1.44269504088896340735992468100189214 +#endif +#ifndef M_LOG10E +# define M_LOG10E 0.434294481903251827651128918916605082 +#endif +#ifndef M_LN2 +# define M_LN2 0.693147180559945309417232121458176568 +#endif +#ifndef M_LN10 +# define M_LN10 2.30258509299404568401799145468436421 +#endif +#ifndef M_PI +# define M_PI 3.14159265358979323846264338327950288 +#endif +#ifndef M_PI_2 +# define M_PI_2 1.57079632679489661923132169163975144 +#endif +#ifndef M_PI_4 +# define M_PI_4 0.785398163397448309615660845819875721 +#endif +#ifndef M_1_PI +# define M_1_PI 0.318309886183790671537767526745028724 +#endif +#ifndef M_2_PI +# define M_2_PI 0.636619772367581343075535053490057448 +#endif +#ifndef M_2_SQRTPI +# define M_2_SQRTPI 1.12837916709551257389615890312154517 +#endif +#ifndef M_SQRT2 +# define M_SQRT2 1.41421356237309504880168872420969808 +#endif +#ifndef M_SQRT1_2 +# define M_SQRT1_2 0.707106781186547524400844362104849039 +#endif + /* C89 math.h: acos asin atan atan2 ceil cos cosh exp fabs floor fmod frexp ldexp @@ -72,7 +112,7 @@ acosh asinh atanh cbrt copysign cosh erf erfc exp2 expm1 fdim fma fmax fmin fpclassify hypot ilogb isfinite isgreater isgreaterequal isinf isless islessequal islessgreater isnan isnormal isunordered - lgamma log1p log2 logb nan nearbyint nextafter nexttoward remainder + lgamma log1p log2 logb lrint nan nearbyint nextafter nexttoward remainder remquo rint round scalbn signbit sinh tanh tgamma trunc * Configure already (5.21.0) scans for: @@ -83,7 +123,10 @@ /* XXX Add ldiv(), lldiv()? It's C99, but from stdlib.h, not math.h */ -/* XXX The truthiness of acosh() is the canary proxy for all of the +/* XXX Beware old gamma() -- one cannot know whether that is the + gamma or the log of gamma, that's why the new tgamma and lgamma. */ + +/* XXX The truthiness of acosh() is the canary for all of the * C99 math. This is very likely wrong, especially in non-UNIX lands * like Win32 and VMS, but also older UNIXes have issues. For Win32 * we later do some undefines for these interfaces. @@ -217,6 +260,7 @@ And the Bessel functions are defined like _this. */ + #ifdef WIN32 # undef c99_exp2 # undef c99_fdim @@ -233,9 +277,18 @@ # undef c99_signbit # undef c99_tgamma # undef c99_trunc + +/* Some APIs exist under Win32 with "underbar" names. */ +# undef c99_hypot +# undef c99_logb +# undef c99_nextafter +# define c99_hypot _hypot +# define c99_logb _logb +# define c99_nextafter _nextafter + #endif -/* The Bessel functions. */ +/* The Bessel functions: BSD, SVID, XPG4, and POSIX. But not C99. */ #ifdef HAS_J0 # if defined(USE_LONG_DOUBLE) && defined(HAS_J0L) # define bessel_j0 j0l @@ -263,12 +316,188 @@ # endif #endif -/* XXX Some of the C99 math functions, if missing, could be rather - * trivially emulated: cbrt, exp2, log2, round, trunc... +/* Emulations for missing math APIs. * * Keep in mind that the point of many of these functions is that * they, if available, are supposed to give more precise/more - * numerically stable results. */ + * numerically stable results. + * + * See e.g. http://www.johndcook.com/math_h.html + */ + +/* XXX cosh sinh tanh */ + +#ifndef c99_acosh +static NV my_acosh(NV x) +{ + return Perl_log(x + Perl_sqrt(x * x - 1)); +} +# define c99_acosh my_acosh +#endif + +#ifndef c99_asinh +static NV my_asinh(NV x) +{ + return Perl_log(x + Perl_sqrt(x * x + 1)); +} +# define c99_asinh my_asinh +#endif + +#ifndef c99_atanh +static NV my_atanh(NV x) +{ + return (Perl_log(1 + x) - Perl_log(1 - x)) / 2; +} +# define c99_atanh my_atanh +#endif + +#ifndef c99_cbrt +static NV my_cbrt(NV x) +{ + static const NV one_third = (NV)1.0/3; + return x >= 0.0 ? Perl_pow(x, one_third) : -Perl_pow(-x, one_third); +} +# define c99_cbrt my_cbrt +#endif + +#ifndef c99_copysign +static NV my_copysign(NV x, NV y) +{ + return y >= 0 ? (x < 0 ? -x : x) : (x < 0 ? x : -x); +} +# define c99_copysign my_copysign +#endif + +#ifndef c99_exp2 +static NV my_exp2(NV x) +{ + return Perl_pow((NV)2.0, x); +} +# define c99_exp2 my_exp2 +#endif + +#ifndef c99_expm1 +static NV my_expm1(NV x) +{ + if (PERL_ABS(x) < 1e-5) + /* Probably not enough for long doubles. */ + return x * (1.0 + x * (0.5 + x / 6.0)); /* Taylor series */ + else + return Perl_exp(x) - 1; +} +# define c99_expm1 my_expm1 +#endif + +#ifndef c99_fdim +static NV my_fdim(NV x, NV y) +{ + return x > y ? x - y : 0; +} +# define c99_fdim my_fdim +#endif + +#ifndef c99_fmax +static NV my_fmax(NV x, NV y) +{ + if (Perl_isnan(x)) { + return Perl_isnan(y) ? NV_NAN : y; + } else if (Perl_isnan(y)) { + return x; + } + return x > y ? x : y; +} +# define c99_fmax my_fmax +#endif + +#ifndef c99_fmin +static NV my_fmin(NV x, NV y) +{ + if (Perl_isnan(x)) { + return Perl_isnan(y) ? NV_NAN : y; + } else if (Perl_isnan(y)) { + return x; + } + return x < y ? x : y; +} +# define c99_fmin my_fmin +#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); + } + return NV_NAN; +} +# define c99_hypot my_hypot +#endif + +#ifndef c99_ilogb +static IV my_ilogb(NV x) +{ + return (IV)(Perl_log(x) * M_LOG2E); +} +# define c99_ilogb my_ilogb +#endif + +#ifndef c99_log1p +static NV my_log1p(NV x) +{ + if (PERL_ABS(x) > 1e-4) + return Perl_log(1.0 + x); + else + /* Probably not enough for long doubles. */ + return x * (1.0 - x * (-x / 2.0 + x / 3.0)); /* Taylor series */ +} +# define c99_log1p my_log1p +#endif + +#ifndef c99_log2 +static NV my_log2(NV x) +{ + return Perl_log(x) * M_LOG2E; +} +# define c99_log2 my_log2 +#endif + +#ifndef c99_round +static NV my_round(NV x) +{ + return (NV)((IV)(x >= 0.0 ? x + 0.5 : x - 0.5)); +} +# define c99_round my_round +#endif + +#ifndef c99_scalbn +# if defined(Perl_ldexpl) && FLT_RADIX == 2 +static NV my_scalbn(NV x) +{ + return Perl_ldexp(x, y); +} +# define c99_scalbn my_scalbn +# endif +#endif + +#ifndef c99_trunc +static NV my_trunc(NV x) +{ + return (NV)((IV)(x)); +} +# define c99_trunc my_trunc +#endif + +#ifndef c99_tgamma +# ifdef c99_lgamma +static NV my_tgamma(NV x) +{ + double l = c99_lgamma(x); + return signgam * Perl_exp(l); /* XXX evil global signgam, need lgamma_r */ +} +# define c99_tgamma my_tgamma +# endif +#endif /* XXX This comment is just to make I_TERMIO and I_SGTTY visible to metaconfig for future extension writers. We don't use them in POSIX. diff --git a/ext/POSIX/t/math.t b/ext/POSIX/t/math.t index e0f5705016..c94984338d 100644 --- a/ext/POSIX/t/math.t +++ b/ext/POSIX/t/math.t @@ -61,11 +61,17 @@ SKIP: { if ($^O =~ /Win32|VMS/) { skip "running in $^O, C99 math support uneven", 28; } - cmp_ok(abs(M_PI - 3.14159265358979), '<', 1e9, "M_PI"); - cmp_ok(abs(asinh(1) - 0.881373587019543), '<', 1e9, "asinh"); - cmp_ok(abs(cbrt(8) - 2), '<', 1e9, "cbrt"); + cmp_ok(abs(M_PI - 3.14159265358979), '<', 1e-9, "M_PI"); + cmp_ok(abs(asinh(1) - 0.881373587019543), '<', 1e-9, "asinh"); + 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"); - cmp_ok(abs(expm1(2) - 6.38905609893065), '<', 1e9, "expm1"); + cmp_ok(abs(expm1(2) - 6.38905609893065), '<', 1e-9, "expm1"); + cmp_ok(abs(expm1(1e-6) - 1.00000050000017e-06), '<', 1e-9, "expm1"); + is(fdim(12, 34), 0, "fdim 12 34"); + is(fdim(34, 12), 22, "fdim 34 12"); + is(fmax(12, 34), 34, "fmax 12 34"); + is(fmin(12, 34), 12, "fmin 12 34"); SKIP: { unless ($Config{d_fpclassify}) { skip "no fpclassify", 4; @@ -75,6 +81,9 @@ SKIP: { is(fpclassify(INFINITY), FP_INFINITE, "fpclassify INFINITY"); is(fpclassify(NAN), FP_NAN, "fpclassify NAN"); } + is(hypot(3, 4), 5, "hypot 3 4"); + is(ilogb(255), 7, "ilogb 255"); + is(ilogb(256), 8, "ilogb 256"); SKIP: { unless ($Config{d_isfinite}) { skip "no isfinite", 1; @@ -93,8 +102,9 @@ SKIP: { } ok(isnan(NAN), "isnan"); } - cmp_ok(abs(log1p(2) - 1.09861228866811), '<', 1e9, "log1p"); - cmp_ok(abs(log2(8) - 3), '<', 1e9, "log2"); + 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"); SKIP: { unless ($Config{d_signbit}) { skip "no signbit", 2; |