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