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
*/
/* 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
/* XXX Regarding C99 math.h, VMS seems to be missing these:
- nan nearbyint round scalbn
+ nan nearbyint round scalbn llrint
*/
#ifdef __VMS
# 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:
/* 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;
+ 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
-static NV my_fpclassify(NV x)
+#ifndef c99_fpclassify
+
+static IV my_fpclassify(NV x)
{
-#if defined(HAS_FPCLASSIFY) && defined(FP_PLUS_INF) /* E.g. HP-UX */
- switch (Perl_fp_class(x)) {
- case FP_PLUS_INF: case FP_MINUS_INF: return FP_INFINITE;
- case FP_SNAN: case FP_QNAN: return FP_NAN;
- case FP_PLUS_NORM: case FP_MINUS_NORM: return FP_NORMAL;
- case FP_PLUS_DENORM: case FP_MINUS_DENORM: return FP_SUBNORMAL;
- case FP_PLUS_ZERO: case FP_MINUS_ZERO: return FP_ZERO;
- default: return -1;
- }
-# define c99_fpclassify my_fpclassify
-#elif (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
-#else
- return -1;
#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
#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
#ifndef c99_scalbn
# if defined(Perl_ldexp) && FLT_RADIX == 2
-static NV my_scalbn(NV x)
+static NV my_scalbn(NV x, int y)
{
return Perl_ldexp(x, y);
}
#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
char* s;
CODE:
#ifdef c99_nan
- RETVAL = c99_nan(s);
+ RETVAL = c99_nan(s ? s : "");
#else
RETVAL = NV_NAN;
not_here("nan");
https://perl5.git.perl.org / perl5.git / blobdiff