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 <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             : 
      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