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

 Generated by: LCOV version 1.10