Python/pymath.c

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