sysdeps/ieee754/ldbl-128/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, write to the Free Software
  61     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA
  62
  63  */
  64
  65 #include "math.h"
  66 #include "math_private.h"
  67
  68 /* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
  69  * 1/sqrt(2) <= x < sqrt(2)
  70  * Theoretical peak relative error = 5.3e-37,
  71  * relative peak error spread = 2.3e-14
  72  */
  73 static const long double P[13] =
  74 {
  75   1.313572404063446165910279910527789794488E4L,
  76   7.771154681358524243729929227226708890930E4L,
  77   2.014652742082537582487669938141683759923E5L,
  78   3.007007295140399532324943111654767187848E5L,
  79   2.854829159639697837788887080758954924001E5L,
  80   1.797628303815655343403735250238293741397E5L,
  81   7.594356839258970405033155585486712125861E4L,
  82   2.128857716871515081352991964243375186031E4L,
  83   3.824952356185897735160588078446136783779E3L,
  84   4.114517881637811823002128927449878962058E2L,
  85   2.321125933898420063925789532045674660756E1L,
  86   4.998469661968096229986658302195402690910E-1L,
  87   1.538612243596254322971797716843006400388E-6L
  88 };
  89 static const long double Q[12] =
  90 {
  91   3.940717212190338497730839731583397586124E4L,
  92   2.626900195321832660448791748036714883242E5L,
  93   7.777690340007566932935753241556479363645E5L,
  94   1.347518538384329112529391120390701166528E6L,
  95   1.514882452993549494932585972882995548426E6L,
  96   1.158019977462989115839826904108208787040E6L,
  97   6.132189329546557743179177159925690841200E5L,
  98   2.248234257620569139969141618556349415120E5L,
  99   5.605842085972455027590989944010492125825E4L,
 100   9.147150349299596453976674231612674085381E3L,
 101   9.104928120962988414618126155557301584078E2L,
 102   4.839208193348159620282142911143429644326E1L
 103 /* 1.000000000000000000000000000000000000000E0L, */
 104 };
 105
 106 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
 107  * where z = 2(x-1)/(x+1)
 108  * 1/sqrt(2) <= x < sqrt(2)
 109  * Theoretical peak relative error = 1.1e-35,
 110  * relative peak error spread 1.1e-9
 111  */
 112 static const long double R[6] =
 113 {
 114   1.418134209872192732479751274970992665513E5L,
 115  -8.977257995689735303686582344659576526998E4L,
 116   2.048819892795278657810231591630928516206E4L,
 117  -2.024301798136027039250415126250455056397E3L,
 118   8.057002716646055371965756206836056074715E1L,
 119  -8.828896441624934385266096344596648080902E-1L
 120 };
 121 static const long double S[6] =
 122 {
 123   1.701761051846631278975701529965589676574E6L,
 124  -1.332535117259762928288745111081235577029E6L,
 125   4.001557694070773974936904547424676279307E5L,
 126  -5.748542087379434595104154610899551484314E4L,
 127   3.998526750980007367835804959888064681098E3L,
 128  -1.186359407982897997337150403816839480438E2L
 129 /* 1.000000000000000000000000000000000000000E0L, */
 130 };
 131
 132 static const long double
 133 /* log10(2) */
 134 L102A = 0.3125L,
 135 L102B = -1.14700043360188047862611052755069732318101185E-2L,
 136 /* log10(e) */
 137 L10EA = 0.5L,
 138 L10EB = -6.570551809674817234887108108339491770560299E-2L,
 139 /* sqrt(2)/2 */
 140 SQRTH = 7.071067811865475244008443621048490392848359E-1L;
 141
 142
 143
 144 /* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 145
 146 static long double
 147 neval (long double x, const long double *p, int n)
 148 {
 149   long double y;
 150
 151   p += n;
 152   y = *p--;
 153   do
 154     {
 155       y = y * x + *p--;
 156     }
 157   while (--n > 0);
 158   return y;
 159 }
 160
 161
 162 /* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 163
 164 static long double
 165 deval (long double x, const long double *p, int n)
 166 {
 167   long double y;
 168
 169   p += n;
 170   y = x + *p--;
 171   do
 172     {
 173       y = y * x + *p--;
 174     }
 175   while (--n > 0);
 176   return y;
 177 }
 178
 179
 180
 181 long double
 182 __ieee754_log10l (long double x)
 183 {
 184   long double z;
 185   long double y;
 186   int e;
 187   int64_t hx, lx;
 188
 189 /* Test for domain */
 190   GET_LDOUBLE_WORDS64 (hx, lx, x);
 191   if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
 192     return (-1.0L / (x - x));
 193   if (hx < 0)
 194     return (x - x) / (x - x);
 195   if (hx >= 0x7fff000000000000LL)
 196     return (x + x);
 197
 198 /* separate mantissa from exponent */
 199
 200 /* Note, frexp is used so that denormal numbers
 201  * will be handled properly.
 202  */
 203   x = __frexpl (x, &e);
 204
 205
 206 /* logarithm using log(x) = z + z**3 P(z)/Q(z),
 207  * where z = 2(x-1)/x+1)
 208  */
 209   if ((e > 2) || (e < -2))
 210     {
 211       if (x < SQRTH)
 212         {                       /* 2( 2x-1 )/( 2x+1 ) */
 213           e -= 1;
 214           z = x - 0.5L;
 215           y = 0.5L * z + 0.5L;
 216         }
 217       else
 218         {                       /*  2 (x-1)/(x+1)   */
 219           z = x - 0.5L;
 220           z -= 0.5L;
 221           y = 0.5L * x + 0.5L;
 222         }
 223       x = z / y;
 224       z = x * x;
 225       y = x * (z * neval (z, R, 5) / deval (z, S, 5));
 226       goto done;
 227     }
 228
 229
 230 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
 231
 232   if (x < SQRTH)
 233     {
 234       e -= 1;
 235       x = 2.0 * x - 1.0L;       /*  2x - 1  */
 236     }
 237   else
 238     {
 239       x = x - 1.0L;
 240     }
 241   z = x * x;
 242   y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
 243   y = y - 0.5 * z;
 244
 245 done:
 246
 247   /* Multiply log of fraction by log10(e)
 248    * and base 2 exponent by log10(2).
 249    */
 250   z = y * L10EB;
 251   z += x * L10EB;
 252   z += e * L102B;
 253   z += y * L10EA;
 254   z += x * L10EA;
 255   z += e * L102A;
 256   return (z);
 257 }
 258 strong_alias (__ieee754_log10l, __log10l_finite)