X-Git-Url: https://perl5.git.perl.org/perl5.git/blobdiff_plain/9d233e1229b9d4643a30bb6859b534d2052d00cd..c36075711e5fbe1f39d9cd2a98188b063f104be6:/ext/POSIX/POSIX.xs

diff --git a/ext/POSIX/POSIX.xs b/ext/POSIX/POSIX.xs
index 5073733..74c0483 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
 
 */
 
@@ -158,12 +158,13 @@
 /* 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
@@ -188,7 +189,7 @@
 #    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
@@ -253,6 +254,23 @@
 #    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
@@ -331,7 +349,7 @@
 
 /* XXX Regarding C99 math.h, VMS seems to be missing these:
 
-  nan nearbyint round scalbn
+  nan nearbyint round scalbn llrint
  */
 
 #ifdef __VMS
@@ -339,6 +357,10 @@
 #    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:
@@ -459,8 +481,36 @@ 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;
+  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)
@@ -474,6 +524,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
@@ -516,83 +568,41 @@ static NV my_fmin(NV x, NV y)
 #  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
@@ -610,6 +620,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
@@ -705,7 +717,7 @@ static NV my_round(NV x)
 
 #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);
 }
@@ -1902,7 +1914,9 @@ acos(x)
 #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
@@ -1964,7 +1978,9 @@ acos(x)
 	    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
@@ -2280,7 +2296,7 @@ nan(s = 0)
 	char*	s;
     CODE:
 #ifdef c99_nan
-	RETVAL = c99_nan(s);
+	RETVAL = c99_nan(s ? s : "");
 #else
 	RETVAL = NV_NAN;
 	not_here("nan");