AppPkg/Applications/Python/Python-2.7.10/Modules/_math.c

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