#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
+#if defined(USE_QUADMATH) && defined(I_QUADMATH)
+
+# undef M_E
+# undef M_LOG2E
+# undef M_LOG10E
+# undef M_LN2
+# undef M_LN10
+# undef M_PI
+# undef M_PI_2
+# undef M_PI_4
+# undef M_1_PI
+# undef M_2_PI
+# undef M_2_SQRTPI
+# undef M_SQRT2
+# undef M_SQRT1_2
+
+# define M_E M_Eq
+# define M_LOG2E M_LOG2Eq
+# define M_LOG10E M_LOG10Eq
+# define M_LN2 M_LN2q
+# define M_LN10 M_LN10q
+# define M_PI M_PIq
+# define M_PI_2 M_PI_2q
+# define M_PI_4 M_PI_4q
+# define M_1_PI M_1_PIq
+# define M_2_PI M_2_PIq
+# define M_2_SQRTPI M_2_SQRTPIq
+# define M_SQRT2 M_SQRT2q
+# define M_SQRT1_2 M_SQRT1_2q
+
+#else
+
+# 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
+
#endif
-#ifndef M_2_PI
-# define M_2_PI 0.636619772367581343075535053490057448
+
+#if !defined(INFINITY) && defined(NV_INF)
+# define INFINITY NV_INF
#endif
-#ifndef M_2_SQRTPI
-# define M_2_SQRTPI 1.12837916709551257389615890312154517
+
+#if !defined(NAN) && defined(NV_NAN)
+# define NAN NV_NAN
#endif
-#ifndef M_SQRT2
-# define M_SQRT2 1.41421356237309504880168872420969808
+
+#if !defined(Inf) && defined(NV_INF)
+# define Inf NV_INF
#endif
-#ifndef M_SQRT1_2
-# define M_SQRT1_2 0.707106781186547524400844362104849039
+
+#if !defined(NaN) && defined(NV_NAN)
+# define NaN NV_NAN
#endif
/* We will have an emulation. */
-#if !defined(HAS_FPCLASSIFY) && !defined(FP_INFINITE)
+#ifndef FP_INFINITE
# define FP_INFINITE 0
# define FP_NAN 1
# define FP_NORMAL 2
# define FP_ZERO 4
#endif
+/* We will have an emulation. */
+#ifndef FE_TONEAREST
+# define FE_TONEAREST 0
+# define FE_TOWARDZERO 1
+# define FE_DOWNWARD 2
+# define FE_UPWARD 3
+#endif
+
/* C89 math.h:
acos asin atan atan2 ceil cos cosh exp fabs floor fmod frexp ldexp
atan2 cos exp log pow sin sqrt
- * Berkeley/SVID extensions:
+ * C99 math.h added:
- j0 j1 jn y0 y1 yn
+ acosh asinh atanh cbrt copysign 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 lrint lround nan nearbyint nextafter nexttoward remainder
+ remquo rint round scalbn signbit tgamma trunc
- * C99 math.h added:
+ See:
+ http://pubs.opengroup.org/onlinepubs/009695399/basedefs/math.h.html
+
+ * Berkeley/SVID extensions:
+
+ j0 j1 jn y0 y1 yn
- 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 lrint nan nearbyint nextafter nexttoward remainder
- remquo rint round scalbn signbit sinh tanh tgamma trunc
+ * 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)
+
+ fegetround fesetround
*/
+/* XXX Constant FP_FAST_FMA (if true, FMA is faster) */
+
/* XXX Add ldiv(), lldiv()? It's C99, but from stdlib.h, not math.h */
/* 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. */
+ * gamma or the log of gamma, that's why the new tgamma and lgamma.
+ * Though also remember tgamma_r and lgamma_r. */
/* 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.
+ * like Win32 and VMS, but also older UNIXes have issues. For Win32,
+ * and other non-fully-C99, we later do some undefines for these interfaces.
*
* But we are very trying very hard to avoid introducing separate Configure
* symbols for all the 40-ish new math symbols. Especially since the set
* of missing functions doesn't seem to follow any patterns. */
#ifdef HAS_ACOSH
-# if defined(USE_LONG_DOUBLE) && defined(HAS_ILOGBL)
+
+/* Certain AIX releases have the C99 math, but not in long double.
+ * The <math.h> has them, e.g. __expl128, but no library has them!
+ *
+ * See the comments in hints/aix.sh about long doubles.
+ *
+ * AIX 5 releases before 5.3 unknown, AIX releases 7 unknown */
+# if defined(_AIX53) || defined(_AIX61)
+# define NO_C99_LONG_DOUBLE_MATH
+# endif
+
+# if defined(USE_QUADMATH) && defined(I_QUADMATH)
+# define c99_acosh acoshq
+# define c99_asinh asinhq
+# define c99_atanh atanhq
+# define c99_cbrt cbrtq
+# define c99_copysign copysignq
+# define c99_erf erfq
+# define c99_erfc erfcq
+/* no exp2q */
+# define c99_expm1 expm1q
+# define c99_fdim fdimq
+# define c99_fma fmaq
+# define c99_fmax fmaxq
+# define c99_fmin fminq
+# define c99_hypot hypotq
+# define c99_ilogb ilogbq
+# define c99_lgamma lgammaq
+# define c99_log1p log1pq
+# define c99_log2 log2q
+/* no logbq */
+/* no llrintq */
+/* no llroundq */
+# define c99_lrint lrintq
+# define c99_lround lroundq
+# define c99_nan nanq
+# define c99_nearbyint nearbyintq
+# define c99_nextafter nextafterq
+/* no nexttowardq */
+# define c99_remainder remainderq
+# define c99_remquo remquoq
+# define c99_rint rintq
+# define c99_round roundq
+# define c99_scalbn scalbnq
+# define c99_signbit signbitq
+# define c99_tgamma tgammal
+# define c99_trunc truncq
+# define bessel_j0 j0q
+# define bessel_j1 j1q
+# define bessel_jn jnq
+# define bessel_y0 y0q
+# define bessel_y1 y1q
+# define bessel_yn ynq
+# elif defined(USE_LONG_DOUBLE) && \
+ !defined(NO_C99_LONG_DOUBLE_MATH) && \
+ defined(HAS_ILOGBL)
/* There's already a symbol for ilogbl, we will use its truthiness
* as the canary for all the *l variants being defined. */
# define c99_acosh acoshl
# define c99_fmin fminl
# define c99_hypot hypotl
# define c99_ilogb ilogbl
-# define c99_lgamma gammal
+# define c99_lgamma lgammal
# define c99_log1p log1pl
# define c99_log2 log2l
# define c99_logb logbl
# else
# define c99_lrint lrintl
# endif
+# if defined(USE_64_BIT_INT) && QUADKIND == QUAD_IS_LONG_LONG
+# define c99_lround llroundl
+# else
+# define c99_lround lroundl
+# endif
# define c99_nan nanl
# define c99_nearbyint nearbyintl
# define c99_nextafter nextafterl
# else
# define c99_lrint lrint
# endif
+# if defined(USE_64_BIT_INT) && QUADKIND == QUAD_IS_LONG_LONG
+# define c99_lround llround
+# else
+# define c99_lround lround
+# endif
# define c99_nan nan
# define c99_nearbyint nearbyint
# define c99_nextafter nextafter
# define c99_trunc trunc
# endif
+# ifndef isunordered
+# ifdef Perl_isnan
+# define isunordered(x, y) (Perl_isnan(x) || Perl_isnan(y))
+# elif defined(HAS_UNORDERED)
+# define isunordered(x, y) unordered(x, y)
+# endif
+# endif
+
+# if !defined(isgreater) && defined(isunordered)
+# define isgreater(x, y) (!isunordered((x), (y)) && (x) > (y))
+# define isgreaterequal(x, y) (!isunordered((x), (y)) && (x) >= (y))
+# define isless(x, y) (!isunordered((x), (y)) && (x) < (y))
+# define islessequal(x, y) (!isunordered((x), (y)) && (x) <= (y))
+# define islessgreater(x, y) (!isunordered((x), (y)) && \
+ ((x) > (y) || (y) > (x)))
+# endif
+
/* Check both the Configure symbol and the macro-ness (like C99 promises). */
# if defined(HAS_FPCLASSIFY) && defined(fpclassify)
# define c99_fpclassify fpclassify
/* If on legacy platforms, and not using gcc, some C99 math interfaces
* might be missing, turn them off so that the emulations hopefully
* kick in. This is admittedly nasty, and fragile, but the alternative
- * is to have Configure scans for all the 40+ interfaces. */
-#ifndef __GNUC__
+ * is to have Configure scans for all the 40+ interfaces.
+ *
+ * For some platforms, also the gcc implementations are missing
+ * certain interfaces.
+ *
+ * In other words: if you have an incomplete (or broken) C99 math interface,
+ * #undef the c99_foo here, and let the emulations kick in. */
-/* HP-UX on PA-RISC is missing certain C99 math functions,
- * but on IA64 (Integrity) these do exist. */
-# if defined(__hpux) && defined(__hppa)
-# undef c99_fma
+#ifdef __GNUC__
+
+/* using gcc */
+
+# if defined(__hpux) && (defined(__hppa) || defined(_PA_RISC))
+# undef c99_nexttoward
+# undef c99_tgamma
+# endif
+
+#else
+
+/* not using gcc */
+
+# if defined(_AIX53) || defined(_AIX61) /* AIX 7 has nexttoward */
# undef c99_nexttoward
-# undef c99_tgamma
+# endif
+
+/* HP-UX on PA-RISC is missing certain C99 math functions,
+ * but on IA64 (Integrity) these do exist, and even on
+ * recent enough HP-UX (cc) releases. */
+# if defined(__hpux) && (defined(__hppa) || defined(_PA_RISC))
+/* lowest known release, could be lower */
+# if defined(__HP_cc) && __HP_cc >= 111120
+# undef c99_fma
+# undef c99_nexttoward
+# undef c99_tgamma
+# else
+# undef c99_exp2
+# undef c99_fdim
+# undef c99_fma
+# undef c99_fmax
+# undef c99_fmin
+# undef c99_fpclassify /* hpux 10.20 has fpclassify but different api */
+# undef c99_lrint
+# undef c99_lround
+# undef c99_nan
+# undef c99_nearbyint
+# undef c99_nexttoward
+# undef c99_remquo
+# undef c99_round
+# undef c99_scalbn
+# undef c99_tgamma
+# undef c99_trunc
+# endif
# endif
# if defined(__irix__)
# undef c99_isinf
# undef c99_isunordered
# undef c99_lrint
+# undef c99_lround
# undef c99_nearbyint
# undef c99_nexttoward
# undef c99_remquo
-# undef c99_rint
# undef c99_round
# undef c99_scalbn
# endif
#endif
+/* XXX Regarding C99 math.h, VMS seems to be missing these:
+
+ lround nan nearbyint round scalbn llrint
+ */
+
+#ifdef __VMS
+# undef c99_lround
+# undef c99_nan
+# undef c99_nearbyint
+# undef c99_round
+# undef c99_scalbn
+/* Have lrint but not llrint. */
+# if defined(USE_64_BIT_INT) && QUADKIND == QUAD_IS_LONG_LONG
+# undef c99_lrint
+# endif
+#endif
+
/* XXX Regarding C99 math.h, Win32 seems to be missing these:
- exp2 fdim fma fmax fmin fpclassify ilogb lgamma log1p log2 lrint
+ erf erfc exp2 fdim fma fmax fmin fpclassify ilogb lgamma log1p log2 lrint
remquo rint signbit tgamma trunc
Win32 does seem to have these:
- acosh asinh atanh cbrt copysign cosh erf erfc expm1 hypot log10 nan
+ acosh asinh atanh cbrt copysign cosh expm1 hypot log10 nan
nearbyint nextafter nexttoward remainder round scalbn
And the Bessel functions are defined like _this.
*/
#ifdef WIN32
+# undef c99_erf
+# undef c99_erfc
# undef c99_exp2
# undef c99_fdim
# undef c99_fma
# undef c99_log1p
# undef c99_log2
# undef c99_lrint
+# undef c99_lround
# undef c99_remquo
# undef c99_rint
# undef c99_signbit
# define c99_logb _logb
# define c99_nextafter _nextafter
+# define bessel_j0 _j0
+# define bessel_j1 _j1
+# define bessel_jn _jn
+# define bessel_y0 _y0
+# define bessel_y1 _y1
+# define bessel_yn _yn
+
+#endif
+
+#ifdef __CYGWIN__
+# undef c99_nexttoward
#endif
/* The Bessel functions: BSD, SVID, XPG4, and POSIX. But not C99. */
-#ifdef HAS_J0
+#if defined(HAS_J0) && !defined(bessel_j0)
# if defined(USE_LONG_DOUBLE) && defined(HAS_J0L)
# define bessel_j0 j0l
# define bessel_j1 j1l
# define bessel_y1 y1l
# define bessel_yn ynl
# else
-# ifdef WIN32
-# define bessel_j0 _j0
-# define bessel_j1 _j1
-# define bessel_jn _jn
-# define bessel_y0 _y0
-# define bessel_y1 _y1
-# define bessel_yn _yn
-# else
-# define bessel_j0 j0
-# define bessel_j1 j1
-# define bessel_jn jn
-# define bessel_y0 y0
-# define bessel_y1 y1
-# define bessel_yn yn
-# endif
+# define bessel_j0 j0
+# define bessel_j1 j1
+# define bessel_jn jn
+# define bessel_y0 y0
+# define bessel_y1 y1
+# define bessel_yn yn
# endif
#endif
* See e.g. http://www.johndcook.com/math_h.html
*/
-/* XXX cosh sinh tanh */
-
#ifndef c99_acosh
static NV my_acosh(NV x)
{
# define c99_copysign my_copysign
#endif
+/* XXX cosh (though c89) */
+
+#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;
+ NV t, y;
+ int sign = x < 0 ? -1 : 1; /* Save the sign. */
+ x = PERL_ABS(x);
+
+ /* Abramowitz and Stegun formula 7.1.26 */
+ t = 1.0 / (1.0 + p * x);
+ y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1) * t * Perl_exp(-x*x);
+
+ return sign * y;
+}
+# define c99_erf my_erf
+#endif
+
+#ifndef c99_erfc
+static NV my_erfc(NV x) {
+ /* This is not necessarily numerically stable, but better than nothing. */
+ return 1.0 - c99_erf(x);
+}
+# define c99_erfc my_erfc
+#endif
+
#ifndef c99_exp2
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.
+ * Taylor series, the first four terms (the last term quartic). */
/* Probably not enough for long doubles. */
- return x * (1.0 + x * (0.5 + x / 6.0)); /* Taylor series */
+ return x * (1.0 + x * (1/2.0 + x * (1/6.0 + x/24.0)));
else
return Perl_exp(x) - 1;
}
# define c99_fmin my_fmin
#endif
-#ifndef HAS_FPCLASSIFY
-static NV my_fpclassify(NV x)
+#ifndef c99_fpclassify
+
+static IV my_fpclassify(NV x)
{
-#if (defined(HAS_FPCLASS) || defined(HAS_FPCLASSL)) && defined(FP_CLASS_SNAN)
- switch (Perl_fp_class(x)) {
- case FP_CLASS_NINF: case FP_CLASS_PINF: return FP_INFINITE;
- case FP_CLASS_SNAN: case FP_CLASS_QNAN: return FP_NAN;
- case FP_CLASS_NNORM: case FP_CLASS_PNORM: return FP_NORMAL;
- case FP_CLASS_NDENORM: case FP_CLASS_PDENORM: return FP_SUBNORMAL;
- case FP_CLASS_NZERO: case FP_CLASS_PZERO: return FP_ZERO;
- default: return -1;
- }
-# define c99_fpclassify my_fpclassify
-#elif (defined(HAS_FPCLASS) || defined(HAS_FP_CLASSL)) && defined(FP_SNAN)
- switch (Perl_fp_class(x)) {
- case FP_NINF: case FP_PINF: return FP_INFINITE;
- case FP_SNAN: case FP_QNAN: return FP_NAN;
- case FP_NNORM: case FP_PNORM: return FP_NORMAL;
- case FP_NDENORM: case FP_PDENORM: return FP_SUBNORMAL;
- case FP_NZERO: case FP_PZERO: return FP_ZERO;
- default: return -1;
- }
-# define c99_fpclassify my_fpclassify
-#elif defined(HAS_FP_CLASS) && defined(FP_POS_INF)
- switch (Perl_fp_class(x)) {
- case FP_NEG_INF: case FP_POS_INF: return FP_INFINITE;
- case FP_SNAN: case FP_QNAN: return FP_NAN;
- case FP_NEG_NORM: case FP_POS_NORM: return FP_NORMAL;
- case FP_NEG_DENORM: case FP_POS_DENORM: return FP_SUBNORMAL;
- case FP_NEG_ZERO: case FP_POS_ZERO: return FP_ZERO;
- default: return -1;
- }
-# define c99_fpclassify my_fpclassify
-#elif defined(HAS_CLASS) && defined(FP_PLUS_INF)
- switch (Perl_fp_class(x)) {
- case FP_MINUS_INF: case FP_PLUS_INF: return FP_INFINITE;
- case FP_SNAN: case FP_QNAN: return FP_NAN;
- case FP_MINUS_NORM: case FP_PLUS_NORM: return FP_NORMAL;
- case FP_MINUS_DENORM: case FP_PLUS_DENORM: return FP_SUBNORMAL;
- case FP_MINUS_ZERO: case FP_PLUS_ZERO: return FP_ZERO;
- default: return -1;
- }
-# define c99_fpclassify my_fpclassify
-#elif defined(HAS_FP_CLASSIFY)
- return Perl_fp_class(x);
-# define c99_fpclassify my_fpclassify
-#elif defined(WIN32)
- int fpclass = _fpclass(x);
+#ifdef Perl_fp_class_inf
if (Perl_fp_class_inf(x)) return FP_INFINITE;
if (Perl_fp_class_nan(x)) return FP_NAN;
if (Perl_fp_class_norm(x)) return FP_NORMAL;
if (Perl_fp_class_denorm(x)) return FP_SUBNORMAL;
if (Perl_fp_class_zero(x)) return FP_ZERO;
- return -1;
# define c99_fpclassify my_fpclassify
#endif
+ return -1;
}
+
#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 * Perl_sqrt(1.0 + t * t);
}
# define c99_hypot my_hypot
#endif
# define c99_ilogb my_ilogb
#endif
+/* XXX lgamma -- non-trivial */
+
#ifndef c99_log1p
static NV my_log1p(NV x)
{
+ /* http://www.johndcook.com/cpp_log_one_plus_x.html -- public domain.
+ * Taylor series, the first four terms (the last term quartic). */
+ if (x < -1.0)
+ return NV_NAN;
+ if (x == -1.0)
+ return -NV_INF;
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 */
+ return x * (1.0 + x * (-1/2.0 + x * (1/3.0 - x/4.0)));
}
# define c99_log1p my_log1p
#endif
# define c99_log2 my_log2
#endif
+/* XXX nextafter */
+
+/* XXX nexttoward */
+
+static int my_fegetround()
+{
+#ifdef HAS_FEGETROUND
+ return fegetround();
+#elif defined(FLT_ROUNDS)
+ return FLT_ROUNDS;
+#elif defined(HAS_FPGETROUND)
+ switch (fpgetround()) {
+ default:
+ case FP_RN: return FE_TONEAREST;
+ case FP_RZ: return FE_TOWARDZERO;
+ case FP_RM: return FE_DOWNWARD;
+ case FE_RP: return FE_UPWARD;
+ }
+#else
+ return -1;
+#endif
+}
+
+/* Toward closest integer. */
+#define MY_ROUND_NEAREST(x) ((NV)((IV)((x) >= 0.0 ? (x) + 0.5 : (x) - 0.5)))
+
+/* Toward zero. */
+#define MY_ROUND_TRUNC(x) ((NV)((IV)(x)))
+
+/* Toward minus infinity. */
+#define MY_ROUND_DOWN(x) ((NV)((IV)((x) >= 0.0 ? (x) : (x) - 0.5)))
+
+/* Toward plus infinity. */
+#define MY_ROUND_UP(x) ((NV)((IV)((x) >= 0.0 ? (x) + 0.5 : (x))))
+
+static NV my_rint(NV x)
+{
+#ifdef FE_TONEAREST
+ switch (my_fegetround()) {
+ default:
+ case FE_TONEAREST: return MY_ROUND_NEAREST(x);
+ case FE_TOWARDZERO: return MY_ROUND_TRUNC(x);
+ case FE_DOWNWARD: return MY_ROUND_DOWN(x);
+ case FE_UPWARD: return MY_ROUND_UP(x);
+ }
+#elif defined(HAS_FPGETROUND)
+ switch (fpgetround()) {
+ default:
+ case FP_RN: return MY_ROUND_NEAREST(x);
+ case FP_RZ: return MY_ROUND_TRUNC(x);
+ case FP_RM: return MY_ROUND_DOWN(x);
+ case FE_RP: return MY_ROUND_UP(x);
+ }
+#else
+ return NV_NAN;
+#endif
+}
+
+/* XXX nearbyint() and rint() are not really identical -- but the difference
+ * is messy: nearbyint is defined NOT to raise FE_INEXACT floating point
+ * exceptions, while rint() is defined to MAYBE raise them. At the moment
+ * Perl is blissfully unaware of such fine detail of floating point. */
+#ifndef c99_nearbyint
+# ifdef FE_TONEAREST
+# define c99_nearbyrint my_rint
+# endif
+#endif
+
+#ifndef c99_lrint
+# ifdef FE_TONEAREST
+static IV my_lrint(NV x)
+{
+ return (IV)my_rint(x);
+}
+# define c99_lrint my_lrint
+# endif
+#endif
+
+#ifndef c99_lround
+static IV my_lround(NV x)
+{
+ return (IV)MY_ROUND_NEAREST(x);
+}
+# define c99_lround my_lround
+#endif
+
+/* XXX remainder */
+
+/* XXX remquo */
+
+#ifndef c99_rint
+# ifdef FE_TONEAREST
+# define c99_rint my_rint
+# endif
+#endif
+
#ifndef c99_round
static NV my_round(NV x)
{
- return (NV)((IV)(x >= 0.0 ? x + 0.5 : x - 0.5));
+ return MY_ROUND_NEAREST(x);
}
# define c99_round my_round
#endif
#ifndef c99_scalbn
-# if defined(Perl_ldexpl) && FLT_RADIX == 2
-static NV my_scalbn(NV x)
+# if defined(Perl_ldexp) && FLT_RADIX == 2
+static NV my_scalbn(NV x, int y)
{
return Perl_ldexp(x, y);
}
# endif
#endif
-#ifndef c99_trunc
-static NV my_trunc(NV x)
-{
- return (NV)((IV)(x));
-}
-# define c99_trunc my_trunc
-#endif
+/* XXX sinh (though c89) */
#ifndef c99_tgamma
# ifdef c99_lgamma
double l = c99_lgamma(x);
return signgam * Perl_exp(l); /* XXX evil global signgam, need lgamma_r */
}
-# define c99_tgamma my_tgamma
+# define c99_tgamma my_tgamma
+/* XXX tgamma without lgamma -- non-trivial */
# endif
#endif
+/* XXX tanh (though c89) */
+
+#ifndef c99_trunc
+static NV my_trunc(NV x)
+{
+ return MY_ROUND_TRUNC(x);
+}
+# define c99_trunc my_trunc
+#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.
(This is really sneaky :-) --AD
y0 = 29
y1 = 30
CODE:
+ RETVAL = NV_NAN;
switch (ix) {
case 0:
- RETVAL = acos(x); /* C89 math */
+ RETVAL = Perl_acos(x); /* C89 math */
break;
case 1:
#ifdef c99_acosh
#endif
break;
case 2:
- RETVAL = asin(x); /* C89 math */
+ RETVAL = Perl_asin(x); /* C89 math */
break;
case 3:
#ifdef c99_asinh
#endif
break;
case 4:
- RETVAL = atan(x); /* C89 math */
+ RETVAL = Perl_atan(x); /* C89 math */
break;
case 5:
#ifdef c99_atanh
#endif
break;
case 7:
- RETVAL = ceil(x); /* C89 math */
+ RETVAL = Perl_ceil(x); /* C89 math */
break;
case 8:
- RETVAL = cosh(x); /* C89 math */
+ RETVAL = Perl_cosh(x); /* C89 math */
break;
case 9:
#ifdef c99_erf
break;
case 10:
#ifdef c99_erfc
- RETVAL = erfc(x);
+ RETVAL = c99_erfc(x);
#else
not_here("erfc");
#endif
#endif
break;
case 13:
- RETVAL = floor(x); /* C89 math */
+ RETVAL = Perl_floor(x); /* C89 math */
break;
case 14:
#ifdef bessel_j0
RETVAL = bessel_j0(x);
#else
- not_here("bessel_j0");
+ not_here("j0");
#endif
break;
case 15:
#ifdef bessel_j1
RETVAL = bessel_j1(x);
#else
- not_here("bessel_j1");
+ not_here("j1");
#endif
break;
case 16:
- /* XXX lgamma_r */
+ /* XXX lgamma_r -- the lgamma accesses a global variable (signgam),
+ * which is evil. Some platforms have lgamma_r, which has
+ * extra parameter instead of the global variable. */
#ifdef c99_lgamma
RETVAL = c99_lgamma(x);
#else
#endif
break;
case 24:
- RETVAL = sinh(x); /* C89 math */
+ RETVAL = Perl_sinh(x); /* C89 math */
break;
case 25:
- RETVAL = tan(x); /* C89 math */
+ RETVAL = Perl_tan(x); /* C89 math */
break;
case 26:
- RETVAL = tanh(x); /* C89 math */
+ RETVAL = Perl_tanh(x); /* C89 math */
break;
case 27:
- /* XXX tgamma_r */
+ /* XXX tgamma_r -- the lgamma accesses a global variable (signgam),
+ * which is evil. Some platforms have tgamma_r, which has
+ * extra parameter instead of the global variable. */
#ifdef c99_tgamma
RETVAL = c99_tgamma(x);
#else
#ifdef bessel_y0
RETVAL = bessel_y0(x);
#else
- not_here("bessel_y0");
+ not_here("y0");
#endif
break;
case 30:
#ifdef bessel_y1
RETVAL = bessel_y1(x);
#else
- not_here("bessel_y1");
+ not_here("y1");
#endif
}
OUTPUT:
fegetround()
CODE:
#ifdef HAS_FEGETROUND
- RETVAL = fegetround();
-#elif defined(FLT_ROUNDS)
- RETVAL = FLT_ROUNDS;
+ RETVAL = my_fegetround();
#else
+ RETVAL = -1;
not_here("fegetround");
#endif
OUTPUT:
CODE:
#ifdef HAS_FEGETROUND /* canary for fesetround */
RETVAL = fesetround(x);
+#elif defined(HAS_FPGETROUND) /* canary for fpsetround */
+ switch (x) {
+ default:
+ case FE_TONEAREST: RETVAL = fpsetround(FP_RN); break;
+ case FE_TOWARDZERO: RETVAL = fpsetround(FP_RZ); break;
+ case FE_DOWNWARD: RETVAL = fpsetround(FP_RM); break;
+ case FE_UPWARD: RETVAL = fpsetround(FP_RP); break;
+ }
#else
+ RETVAL = -1;
not_here("fesetround");
#endif
OUTPUT:
isinf = 3
isnan = 4
isnormal = 5
- signbit = 6
+ lrint = 6
+ lround = 7
+ signbit = 8
CODE:
+ RETVAL = -1;
switch (ix) {
case 0:
#ifdef c99_fpclassify
#endif
break;
case 6:
+#ifdef c99_lrint
+ RETVAL = c99_lrint(x);
+#else
+ not_here("lrint");
+#endif
+ break;
+ case 7:
+#ifdef c99_lround
+ RETVAL = c99_lround(x);
+#else
+ not_here("lround");
+#endif
+ break;
+ case 8:
default:
#ifdef Perl_signbit
RETVAL = Perl_signbit(x);
nexttoward = 13
remainder = 14
CODE:
+ RETVAL = NV_NAN;
switch (ix) {
case 0:
#ifdef c99_copysign
#endif
break;
case 4:
- RETVAL = fmod(x, y); /* C89 math */
+ RETVAL = Perl_fmod(x, y); /* C89 math */
break;
case 5:
#ifdef c99_hypot
#ifdef c99_scalbn
RETVAL = c99_scalbn(x, y);
#else
+ RETVAL = NV_NAN;
not_here("scalbn");
#endif
OUTPUT:
#ifdef c99_fma
RETVAL = c99_fma(x, y, z);
#else
+ RETVAL = NV_NAN;
not_here("fma");
#endif
OUTPUT:
char* s;
CODE:
#ifdef c99_nan
- RETVAL = c99_nan(s);
+ RETVAL = c99_nan(s ? s : "");
#else
+ RETVAL = NV_NAN;
+# ifndef NV_NAN
not_here("nan");
+# endif
#endif
OUTPUT:
RETVAL
ALIAS:
yn = 1
CODE:
+ RETVAL = NV_NAN;
switch (ix) {
case 0:
#ifdef bessel_jn
https://perl5.git.perl.org / perl5.git / blobdiff