sysdeps/ieee754/ldbl-128ibm/e_log10l.c

   1 /*                                                      log10l.c
   2  *
   3  *      Common logarithm, 128-bit long double precision
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * long double x, y, log10l();
  10  *
  11  * y = log10l( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns the base 10 logarithm of x.
  18  *
  19  * The argument is separated into its exponent and fractional
  20  * parts.  If the exponent is between -1 and +1, the logarithm
  21  * of the fraction is approximated by
  22  *
  23  *     log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
  24  *
  25  * Otherwise, setting  z = 2(x-1)/x+1),
  26  *
  27  *     log(x) = z + z^3 P(z)/Q(z).
  28  *
  29  *
  30  *
  31  * ACCURACY:
  32  *
  33  *                      Relative error:
  34  * arithmetic   domain     # trials      peak         rms
  35  *    IEEE      0.5, 2.0     30000      2.3e-34     4.9e-35
  36  *    IEEE     exp(+-10000)  30000      1.0e-34     4.1e-35
  37  *
  38  * In the tests over the interval exp(+-10000), the logarithms
  39  * of the random arguments were uniformly distributed over
  40  * [-10000, +10000].
  41  *
  42  */
  43
  44 /*
  45    Cephes Math Library Release 2.2:  January, 1991
  46    Copyright 1984, 1991 by Stephen L. Moshier
  47    Adapted for glibc November, 2001
  48
  49     This library is free software; you can redistribute it and/or
  50     modify it under the terms of the GNU Lesser General Public
  51     License as published by the Free Software Foundation; either
  52     version 2.1 of the License, or (at your option) any later version.
  53
  54     This library is distributed in the hope that it will be useful,
  55     but WITHOUT ANY WARRANTY; without even the implied warranty of
  56     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  57     Lesser General Public License for more details.
  58
  59     You should have received a copy of the GNU Lesser General Public
  60     License along with this library; if not, see <http://www.gnu.org/licenses/>.
  61  */
  62
  63 #include "math.h"
  64 #include "math_private.h"
  65
  66 /* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
  67  * 1/sqrt(2) <= x < sqrt(2)
  68  * Theoretical peak relative error = 5.3e-37,
  69  * relative peak error spread = 2.3e-14
  70  */
  71 static const long double P[13] =
  72 {
  73   1.313572404063446165910279910527789794488E4L,
  74   7.771154681358524243729929227226708890930E4L,
  75   2.014652742082537582487669938141683759923E5L,
  76   3.007007295140399532324943111654767187848E5L,
  77   2.854829159639697837788887080758954924001E5L,
  78   1.797628303815655343403735250238293741397E5L,
  79   7.594356839258970405033155585486712125861E4L,
  80   2.128857716871515081352991964243375186031E4L,
  81   3.824952356185897735160588078446136783779E3L,
  82   4.114517881637811823002128927449878962058E2L,
  83   2.321125933898420063925789532045674660756E1L,
  84   4.998469661968096229986658302195402690910E-1L,
  85   1.538612243596254322971797716843006400388E-6L
  86 };
  87 static const long double Q[12] =
  88 {
  89   3.940717212190338497730839731583397586124E4L,
  90   2.626900195321832660448791748036714883242E5L,
  91   7.777690340007566932935753241556479363645E5L,
  92   1.347518538384329112529391120390701166528E6L,
  93   1.514882452993549494932585972882995548426E6L,
  94   1.158019977462989115839826904108208787040E6L,
  95   6.132189329546557743179177159925690841200E5L,
  96   2.248234257620569139969141618556349415120E5L,
  97   5.605842085972455027590989944010492125825E4L,
  98   9.147150349299596453976674231612674085381E3L,
  99   9.104928120962988414618126155557301584078E2L,
 100   4.839208193348159620282142911143429644326E1L
 101 /* 1.000000000000000000000000000000000000000E0L, */
 102 };
 103
 104 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
 105  * where z = 2(x-1)/(x+1)
 106  * 1/sqrt(2) <= x < sqrt(2)
 107  * Theoretical peak relative error = 1.1e-35,
 108  * relative peak error spread 1.1e-9
 109  */
 110 static const long double R[6] =
 111 {
 112   1.418134209872192732479751274970992665513E5L,
 113  -8.977257995689735303686582344659576526998E4L,
 114   2.048819892795278657810231591630928516206E4L,
 115  -2.024301798136027039250415126250455056397E3L,
 116   8.057002716646055371965756206836056074715E1L,
 117  -8.828896441624934385266096344596648080902E-1L
 118 };
 119 static const long double S[6] =
 120 {
 121   1.701761051846631278975701529965589676574E6L,
 122  -1.332535117259762928288745111081235577029E6L,
 123   4.001557694070773974936904547424676279307E5L,
 124  -5.748542087379434595104154610899551484314E4L,
 125   3.998526750980007367835804959888064681098E3L,
 126  -1.186359407982897997337150403816839480438E2L
 127 /* 1.000000000000000000000000000000000000000E0L, */
 128 };
 129
 130 static const long double
 131 /* log10(2) */
 132 L102A = 0.3125L,
 133 L102B = -1.14700043360188047862611052755069732318101185E-2L,
 134 /* log10(e) */
 135 L10EA = 0.5L,
 136 L10EB = -6.570551809674817234887108108339491770560299E-2L,
 137 /* sqrt(2)/2 */
 138 SQRTH = 7.071067811865475244008443621048490392848359E-1L;
 139
 140
 141
 142 /* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 143
 144 static long double
 145 neval (long double x, const long double *p, int n)
 146 {
 147   long double y;
 148
 149   p += n;
 150   y = *p--;
 151   do
 152     {
 153       y = y * x + *p--;
 154     }
 155   while (--n > 0);
 156   return y;
 157 }
 158
 159
 160 /* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 161
 162 static long double
 163 deval (long double x, const long double *p, int n)
 164 {
 165   long double y;
 166
 167   p += n;
 168   y = x + *p--;
 169   do
 170     {
 171       y = y * x + *p--;
 172     }
 173   while (--n > 0);
 174   return y;
 175 }
 176
 177
 178
 179 long double
 180 __ieee754_log10l (long double x)
 181 {
 182   long double z;
 183   long double y;
 184   int e;
 185   int64_t hx, lx;
 186
 187 /* Test for domain */
 188   GET_LDOUBLE_WORDS64 (hx, lx, x);
 189   if (((hx & 0x7fffffffffffffffLL) | (lx & 0x7fffffffffffffffLL)) == 0)
 190     return (-1.0L / (x - x));
 191   if (hx < 0)
 192     return (x - x) / (x - x);
 193   if (hx >= 0x7ff0000000000000LL)
 194     return (x + x);
 195
 196 /* separate mantissa from exponent */
 197
 198 /* Note, frexp is used so that denormal numbers
 199  * will be handled properly.
 200  */
 201   x = __frexpl (x, &e);
 202
 203
 204 /* logarithm using log(x) = z + z**3 P(z)/Q(z),
 205  * where z = 2(x-1)/x+1)
 206  */
 207   if ((e > 2) || (e < -2))
 208     {
 209       if (x < SQRTH)
 210         {                       /* 2( 2x-1 )/( 2x+1 ) */
 211           e -= 1;
 212           z = x - 0.5L;
 213           y = 0.5L * z + 0.5L;
 214         }
 215       else
 216         {                       /*  2 (x-1)/(x+1)   */
 217           z = x - 0.5L;
 218           z -= 0.5L;
 219           y = 0.5L * x + 0.5L;
 220         }
 221       x = z / y;
 222       z = x * x;
 223       y = x * (z * neval (z, R, 5) / deval (z, S, 5));
 224       goto done;
 225     }
 226
 227
 228 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
 229
 230   if (x < SQRTH)
 231     {
 232       e -= 1;
 233       x = 2.0 * x - 1.0L;       /*  2x - 1  */
 234     }
 235   else
 236     {
 237       x = x - 1.0L;
 238     }
 239   z = x * x;
 240   y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
 241   y = y - 0.5 * z;
 242
 243 done:
 244
 245   /* Multiply log of fraction by log10(e)
 246    * and base 2 exponent by log10(2).
 247    */
 248   z = y * L10EB;
 249   z += x * L10EB;
 250   z += e * L102B;
 251   z += y * L10EA;
 252   z += x * L10EA;
 253   z += e * L102A;
 254   return (z);
 255 }
 256 strong_alias (__ieee754_log10l, __log10l_finite)