Python/pymath.c

   1 #include "Python.h"
   2
   3 #ifdef X87_DOUBLE_ROUNDING
   4 /* On x86 platforms using an x87 FPU, this function is called from the
   5    Py_FORCE_DOUBLE macro (defined in pymath.h) to force a floating-point
   6    number out of an 80-bit x87 FPU register and into a 64-bit memory location,
   7    thus rounding from extended precision to double precision. */
   8 double _Py_force_double(double x)
   9 {
  10         volatile double y;
  11         y = x;
  12         return y;
  13 }
  14 #endif
  15
  16 #ifndef HAVE_HYPOT
  17 double hypot(double x, double y)
  18 {
  19         double yx;
  20
  21         x = fabs(x);
  22         y = fabs(y);
  23         if (x < y) {
  24                 double temp = x;
  25                 x = y;
  26                 y = temp;
  27         }
  28         if (x == 0.)
  29                 return 0.;
  30         else {
  31                 yx = y/x;
  32                 return x*sqrt(1.+yx*yx);
  33         }
  34 }
  35 #endif /* HAVE_HYPOT */
  36
  37 #ifndef HAVE_COPYSIGN
  38 double
  39 copysign(double x, double y)
  40 {
  41         /* use atan2 to distinguish -0. from 0. */
  42         if (y > 0. || (y == 0. && atan2(y, -1.) > 0.)) {
  43                 return fabs(x);
  44         } else {
  45                 return -fabs(x);
  46         }
  47 }
  48 #endif /* HAVE_COPYSIGN */
  49
  50 #ifndef HAVE_ROUND
  51 double
  52 round(double x)
  53 {
  54     double absx, y;
  55     absx = fabs(x);
  56     y = floor(absx);
  57     if (absx - y >= 0.5)
  58         y += 1.0;
  59     return copysign(y, x);
  60 }
  61 #endif /* HAVE_ROUND */
  62
  63 #ifndef HAVE_LOG1P
  64 #include <float.h>
  65
  66 double
  67 log1p(double x)
  68 {
  69         /* For x small, we use the following approach.  Let y be the nearest
  70            float to 1+x, then
  71
  72              1+x = y * (1 - (y-1-x)/y)
  73
  74            so log(1+x) = log(y) + log(1-(y-1-x)/y).  Since (y-1-x)/y is tiny,
  75            the second term is well approximated by (y-1-x)/y.  If abs(x) >=
  76            DBL_EPSILON/2 or the rounding-mode is some form of round-to-nearest
  77            then y-1-x will be exactly representable, and is computed exactly
  78            by (y-1)-x.
  79
  80            If abs(x) < DBL_EPSILON/2 and the rounding mode is not known to be
  81            round-to-nearest then this method is slightly dangerous: 1+x could
  82            be rounded up to 1+DBL_EPSILON instead of down to 1, and in that
  83            case y-1-x will not be exactly representable any more and the
  84            result can be off by many ulps.  But this is easily fixed: for a
  85            floating-point number |x| < DBL_EPSILON/2., the closest
  86            floating-point number to log(1+x) is exactly x.
  87         */
  88
  89         double y;
  90         if (fabs(x) < DBL_EPSILON/2.) {
  91                 return x;
  92         } else if (-0.5 <= x && x <= 1.) {
  93                 /* WARNING: it's possible than an overeager compiler
  94                    will incorrectly optimize the following two lines
  95                    to the equivalent of "return log(1.+x)". If this
  96                    happens, then results from log1p will be inaccurate
  97                    for small x. */
  98                 y = 1.+x;
  99                 return log(y)-((y-1.)-x)/y;
 100         } else {
 101                 /* NaNs and infinities should end up here */
 102                 return log(1.+x);
 103         }
 104 }
 105 #endif /* HAVE_LOG1P */
 106
 107 /*
 108  * ====================================================
 109  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 110  *
 111  * Developed at SunPro, a Sun Microsystems, Inc. business.
 112  * Permission to use, copy, modify, and distribute this
 113  * software is freely granted, provided that this notice
 114  * is preserved.
 115  * ====================================================
 116  */
 117
 118 static const double ln2 = 6.93147180559945286227E-01;
 119 static const double two_pow_m28 = 3.7252902984619141E-09; /* 2**-28 */
 120 static const double two_pow_p28 = 268435456.0; /* 2**28 */
 121 static const double zero = 0.0;
 122
 123 /* asinh(x)
 124  * Method :
 125  *      Based on
 126  *              asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
 127  *      we have
 128  *      asinh(x) := x  if  1+x*x=1,
 129  *               := sign(x)*(log(x)+ln2)) for large |x|, else
 130  *               := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
 131  *               := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
 132  */
 133
 134 #ifndef HAVE_ASINH
 135 double
 136 asinh(double x)
 137 {
 138         double w;
 139         double absx = fabs(x);
 140
 141         if (Py_IS_NAN(x) || Py_IS_INFINITY(x)) {
 142                 return x+x;
 143         }
 144         if (absx < two_pow_m28) {       /* |x| < 2**-28 */
 145                 return x;       /* return x inexact except 0 */
 146         }
 147         if (absx > two_pow_p28) {       /* |x| > 2**28 */
 148                 w = log(absx)+ln2;
 149         }
 150         else if (absx > 2.0) {          /* 2 < |x| < 2**28 */
 151                 w = log(2.0*absx + 1.0 / (sqrt(x*x + 1.0) + absx));
 152         }
 153         else {                          /* 2**-28 <= |x| < 2= */
 154                 double t = x*x;
 155                 w = log1p(absx + t / (1.0 + sqrt(1.0 + t)));
 156         }
 157         return copysign(w, x);
 158
 159 }
 160 #endif /* HAVE_ASINH */
 161
 162 /* acosh(x)
 163  * Method :
 164  *      Based on
 165  *            acosh(x) = log [ x + sqrt(x*x-1) ]
 166  *      we have
 167  *            acosh(x) := log(x)+ln2, if x is large; else
 168  *            acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
 169  *            acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
 170  *
 171  * Special cases:
 172  *      acosh(x) is NaN with signal if x<1.
 173  *      acosh(NaN) is NaN without signal.
 174  */
 175
 176 #ifndef HAVE_ACOSH
 177 double
 178 acosh(double x)
 179 {
 180         if (Py_IS_NAN(x)) {
 181                 return x+x;
 182         }
 183         if (x < 1.) {                   /* x < 1;  return a signaling NaN */
 184                 errno = EDOM;
 185 #ifdef Py_NAN
 186                 return Py_NAN;
 187 #else
 188                 return (x-x)/(x-x);
 189 #endif
 190         }
 191         else if (x >= two_pow_p28) {    /* x > 2**28 */
 192                 if (Py_IS_INFINITY(x)) {
 193                         return x+x;
 194                 } else {
 195                         return log(x)+ln2;      /* acosh(huge)=log(2x) */
 196                 }
 197         }
 198         else if (x == 1.) {
 199                 return 0.0;                     /* acosh(1) = 0 */
 200         }
 201         else if (x > 2.) {                      /* 2 < x < 2**28 */
 202                 double t = x*x;
 203                 return log(2.0*x - 1.0 / (x + sqrt(t - 1.0)));
 204         }
 205         else {                          /* 1 < x <= 2 */
 206                 double t = x - 1.0;
 207                 return log1p(t + sqrt(2.0*t + t*t));
 208         }
 209 }
 210 #endif /* HAVE_ACOSH */
 211
 212 /* atanh(x)
 213  * Method :
 214  *    1.Reduced x to positive by atanh(-x) = -atanh(x)
 215  *    2.For x>=0.5
 216  *                1           2x                          x
 217  *      atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
 218  *                2          1 - x                    1 - x
 219  *
 220  *      For x<0.5
 221  *      atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
 222  *
 223  * Special cases:
 224  *      atanh(x) is NaN if |x| >= 1 with signal;
 225  *      atanh(NaN) is that NaN with no signal;
 226  *
 227  */
 228
 229 #ifndef HAVE_ATANH
 230 double
 231 atanh(double x)
 232 {
 233         double absx;
 234         double t;
 235
 236         if (Py_IS_NAN(x)) {
 237                 return x+x;
 238         }
 239         absx = fabs(x);
 240         if (absx >= 1.) {               /* |x| >= 1 */
 241                 errno = EDOM;
 242 #ifdef Py_NAN
 243                 return Py_NAN;
 244 #else
 245                 return x/zero;
 246 #endif
 247         }
 248         if (absx < two_pow_m28) {       /* |x| < 2**-28 */
 249                 return x;
 250         }
 251         if (absx < 0.5) {               /* |x| < 0.5 */
 252                 t = absx+absx;
 253                 t = 0.5 * log1p(t + t*absx / (1.0 - absx));
 254         }
 255         else {                          /* 0.5 <= |x| <= 1.0 */
 256                 t = 0.5 * log1p((absx + absx) / (1.0 - absx));
 257         }
 258         return copysign(t, x);
 259 }
 260 #endif /* HAVE_ATANH */