#ifdef I_FLOAT
#include <float.h>
#endif
+#ifdef I_FENV
+#include <fenv.h>
+#endif
#ifdef I_LIMITS
#include <limits.h>
#endif
# define M_SQRT1_2 0.707106781186547524400844362104849039
#endif
+#if !defined(INFINITY) && defined(NV_INF)
+# define INFINITY NV_INF
+#endif
+
+#if !defined(NAN) && defined(NV_NAN)
+# define NAN NV_NAN
+#endif
+
+#if !defined(Inf) && defined(NV_INF)
+# define Inf NV_INF
+#endif
+
+#if !defined(NaN) && defined(NV_NAN)
+# define NaN NV_NAN
+#endif
+
+/* We will have an emulation. */
+#ifndef FP_INFINITE
+# define FP_INFINITE 0
+# define FP_NAN 1
+# define FP_NORMAL 2
+# define FP_SUBNORMAL 3
+# define FP_ZERO 4
+#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:
+
+ 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 nan nearbyint nextafter nexttoward remainder
+ remquo rint round scalbn signbit tgamma trunc
- j0 j1 jn y0 y1 yn
+ * Berkeley/SVID extensions:
- * C99 math.h added:
+ 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
# 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
# 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
# endif
#endif
+/* 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.
+ *
+ * In other words: if you have an incomplete (or broken) C99 math interface,
+ * #undef the c99_foo here, and let the emulations kick in. */
+
+#ifndef __GNUC__
+
/* 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
-# undef c99_nexttoward
-# undef c99_tgamma
+# if defined(__hpux) && defined(__hppa)
+# undef c99_exp2
+# undef c99_fdim
+# undef c99_fma
+# undef c99_fmax
+# undef c99_fmin
+# undef c99_fpclassify
+# undef c99_lrint
+# 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
+
+# if defined(__irix__)
+# undef c99_ilogb
+# undef c99_exp2
+# endif
+
+# if defined(__osf__) /* Tru64 */
+# undef c99_fdim
+# undef c99_fma
+# undef c99_fmax
+# undef c99_fmin
+# undef c99_fpclassify
+# undef c99_isfinite
+# undef c99_isinf
+# undef c99_isunordered
+# undef c99_lrint
+# 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:
+
+ nan nearbyint round scalbn llrint
+ */
+
+#ifdef __VMS
+# 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:
# 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
/* The Bessel functions: BSD, SVID, XPG4, and POSIX. But not C99. */
# 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 * 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.
+ * Also including the cubic term. */
/* Probably not enough for long doubles. */
return x * (1.0 + x * (0.5 + x / 6.0)); /* Taylor series */
else
# define c99_fmin my_fmin
#endif
+#ifndef c99_fpclassify
+
+static IV my_fpclassify(NV 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;
+# 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 * 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.
+ * Including also quadratic term. */
if (PERL_ABS(x) > 1e-4)
return Perl_log(1.0 + x);
else
# 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;
+ /* XXX emulate using fpgetround() (HAS_FPGETROUND):
+ * FP_RN to nearest, FP_RM down, FP_RP, up, FP_RZ truncate */
+#else
+ return -1;
+#endif
+}
+
+static NV my_rint(NV x)
+{
+#ifdef FE_TONEAREST
+ switch (my_fegetround()) {
+ default:
+ case FE_TONEAREST:
+ return (NV)((IV)(x >= 0.0 ? x + 0.5 : x - 0.5)); /* like round() */
+ case FE_TOWARDZERO:
+ return (NV)((IV)(x)); /* like trunc() */
+ case FE_DOWNWARD:
+ return (NV)((IV)(x >= 0.0 ? x : x - 0.5));
+ case FE_UPWARD:
+ return (NV)((IV)(x >= 0.0 ? x + 0.5 : x));
+ }
+#else
+ /* XXX emulate using fpsetround() (HAS_FPGETROUND):
+ * FP_RN to nearest, FP_RM down, FP_RP, up, FP_RZ truncate */
+ 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 lrint(NV x)
+{
+ return (IV)my_rint(x);
+}
+# define c99_lrint my_lrint
+# endif
+#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)
{
#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 (NV)((IV)(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 */
#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
RETVAL = 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:
RETVAL
IV
+fegetround()
+ CODE:
+#ifdef HAS_FEGETROUND
+ RETVAL = my_fegetround();
+#else
+ RETVAL = -1;
+ not_here("fegetround");
+#endif
+ OUTPUT:
+ RETVAL
+
+IV
+fesetround(x)
+ IV x
+ CODE:
+#ifdef HAS_FEGETROUND /* canary for fesetround */
+ RETVAL = fesetround(x);
+#else
+ RETVAL = -1;
+ not_here("fesetround");
+#endif
+ OUTPUT:
+ RETVAL
+
+IV
fpclassify(x)
NV x
ALIAS:
isinf = 3
isnan = 4
isnormal = 5
- signbit = 6
+ lrint = 6
+ signbit = 7
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:
default:
#ifdef Perl_signbit
RETVAL = Perl_signbit(x);
nexttoward = 13
remainder = 14
CODE:
+ RETVAL = NV_NAN;
switch (ix) {
case 0:
#ifdef c99_copysign
#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;
not_here("nan");
#endif
OUTPUT:
ALIAS:
yn = 1
CODE:
+ RETVAL = NV_NAN;
switch (ix) {
case 0:
#ifdef bessel_jn
https://perl5.git.perl.org / perl5.git / blobdiff