nsz Git - libm/blob - src/math/logl.c

   1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_logl.c */
   2 /*
   3  * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
   4  *
   5  * Permission to use, copy, modify, and distribute this software for any
   6  * purpose with or without fee is hereby granted, provided that the above
   7  * copyright notice and this permission notice appear in all copies.
   8  *
   9  * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
  10  * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
  11  * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
  12  * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  13  * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
  14  * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
  15  * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  16  */
  17 /*
  18  *      Natural logarithm, long double precision
  19  *
  20  *
  21  * SYNOPSIS:
  22  *
  23  * long double x, y, logl();
  24  *
  25  * y = logl( x );
  26  *
  27  *
  28  * DESCRIPTION:
  29  *
  30  * Returns the base e (2.718...) logarithm of x.
  31  *
  32  * The argument is separated into its exponent and fractional
  33  * parts.  If the exponent is between -1 and +1, the logarithm
  34  * of the fraction is approximated by
  35  *
  36  *     log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
  37  *
  38  * Otherwise, setting  z = 2(x-1)/x+1),
  39  *
  40  *     log(x) = z + z**3 P(z)/Q(z).
  41  *
  42  *
  43  * ACCURACY:
  44  *
  45  *                      Relative error:
  46  * arithmetic   domain     # trials      peak         rms
  47  *    IEEE      0.5, 2.0    150000      8.71e-20    2.75e-20
  48  *    IEEE     exp(+-10000) 100000      5.39e-20    2.34e-20
  49  *
  50  * In the tests over the interval exp(+-10000), the logarithms
  51  * of the random arguments were uniformly distributed over
  52  * [-10000, +10000].
  53  *
  54  * ERROR MESSAGES:
  55  *
  56  * log singularity:  x = 0; returns -INFINITY
  57  * log domain:       x < 0; returns NAN
  58  */
  59
  60 #include "libm.h"
  61
  62 #if LD64
  63 long double logl(long double x)
  64 {
  65         return log(x);
  66 }
  67 #elif LD80
  68 /* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
  69  * 1/sqrt(2) <= x < sqrt(2)
  70  * Theoretical peak relative error = 2.32e-20
  71  */
  72 static long double P[] = {
  73  4.5270000862445199635215E-5L,
  74  4.9854102823193375972212E-1L,
  75  6.5787325942061044846969E0L,
  76  2.9911919328553073277375E1L,
  77  6.0949667980987787057556E1L,
  78  5.7112963590585538103336E1L,
  79  2.0039553499201281259648E1L,
  80 };
  81 static long double Q[] = {
  82 /* 1.0000000000000000000000E0,*/
  83  1.5062909083469192043167E1L,
  84  8.3047565967967209469434E1L,
  85  2.2176239823732856465394E2L,
  86  3.0909872225312059774938E2L,
  87  2.1642788614495947685003E2L,
  88  6.0118660497603843919306E1L,
  89 };
  90
  91 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
  92  * where z = 2(x-1)/(x+1)
  93  * 1/sqrt(2) <= x < sqrt(2)
  94  * Theoretical peak relative error = 6.16e-22
  95  */
  96 static long double R[4] = {
  97  1.9757429581415468984296E-3L,
  98 -7.1990767473014147232598E-1L,
  99  1.0777257190312272158094E1L,
 100 -3.5717684488096787370998E1L,
 101 };
 102 static long double S[4] = {
 103 /* 1.00000000000000000000E0L,*/
 104 -2.6201045551331104417768E1L,
 105  1.9361891836232102174846E2L,
 106 -4.2861221385716144629696E2L,
 107 };
 108 static const long double C1 = 6.9314575195312500000000E-1L;
 109 static const long double C2 = 1.4286068203094172321215E-6L;
 110
 111 #define SQRTH 0.70710678118654752440L
 112
 113 long double logl(long double x)
 114 {
 115         long double y, z;
 116         int e;
 117
 118         if (isnan(x))
 119                 return x;
 120         if (x == INFINITY)
 121                 return x;
 122         if (x <= 0.0L) {
 123                 if (x == 0.0L)
 124                         return -INFINITY;
 125                 return NAN;
 126         }
 127
 128         /* separate mantissa from exponent */
 129         /* Note, frexp is used so that denormal numbers
 130          * will be handled properly.
 131          */
 132         x = frexpl(x, &e);
 133
 134         /* logarithm using log(x) = z + z**3 P(z)/Q(z),
 135          * where z = 2(x-1)/x+1)
 136          */
 137         if (e > 2 || e < -2) {
 138                 if (x < SQRTH) {  /* 2(2x-1)/(2x+1) */
 139                         e -= 1;
 140                         z = x - 0.5L;
 141                         y = 0.5L * z + 0.5L;
 142                 } else {  /*  2 (x-1)/(x+1)   */
 143                         z = x - 0.5L;
 144                         z -= 0.5L;
 145                         y = 0.5L * x  + 0.5L;
 146                 }
 147                 x = z / y;
 148                 z = x*x;
 149                 z = x * (z * __polevll(z, R, 3) / __p1evll(z, S, 3));
 150                 z = z + e * C2;
 151                 z = z + x;
 152                 z = z + e * C1;
 153                 return z;
 154         }
 155
 156         /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
 157         if (x < SQRTH) {
 158                 e -= 1;
 159                 x = ldexpl(x, 1) - 1.0L; /*  2x - 1  */
 160         } else {
 161                 x = x - 1.0L;
 162         }
 163         z = x*x;
 164         y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 6));
 165         y = y + e * C2;
 166         z = y - ldexpl(z, -1);   /*  y - 0.5 * z  */
 167         /* Note, the sum of above terms does not exceed x/4,
 168          * so it contributes at most about 1/4 lsb to the error.
 169          */
 170         z = z + x;
 171         z = z + e * C1; /* This sum has an error of 1/2 lsb. */
 172         return z;
 173 }
 174 #endif