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 (x)
 183      long double x;
 184 {
 185   long double z;
 186   long double y;
 187   int e;
 188   int64_t hx, lx;
 189
 190 /* Test for domain */
 191   GET_LDOUBLE_WORDS64 (hx, lx, x);
 192   if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
 193     return (-1.0L / (x - x));
 194   if (hx < 0)
 195     return (x - x) / (x - x);
 196   if (hx >= 0x7fff000000000000LL)
 197     return (x + x);
 198
 199 /* separate mantissa from exponent */
 200
 201 /* Note, frexp is used so that denormal numbers
 202  * will be handled properly.
 203  */
 204   x = __frexpl (x, &e);
 205
 206
 207 /* logarithm using log(x) = z + z**3 P(z)/Q(z),
 208  * where z = 2(x-1)/x+1)
 209  */
 210   if ((e > 2) || (e < -2))
 211     {
 212       if (x < SQRTH)
 213         {                       /* 2( 2x-1 )/( 2x+1 ) */
 214           e -= 1;
 215           z = x - 0.5L;
 216           y = 0.5L * z + 0.5L;
 217         }
 218       else
 219         {                       /*  2 (x-1)/(x+1)   */
 220           z = x - 0.5L;
 221           z -= 0.5L;
 222           y = 0.5L * x + 0.5L;
 223         }
 224       x = z / y;
 225       z = x * x;
 226       y = x * (z * neval (z, R, 5) / deval (z, S, 5));
 227       goto done;
 228     }
 229
 230
 231 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
 232
 233   if (x < SQRTH)
 234     {
 235       e -= 1;
 236       x = 2.0 * x - 1.0L;       /*  2x - 1  */
 237     }
 238   else
 239     {
 240       x = x - 1.0L;
 241     }
 242   z = x * x;
 243   y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
 244   y = y - 0.5 * z;
 245
 246 done:
 247
 248   /* Multiply log of fraction by log10(e)
 249    * and base 2 exponent by log10(2).
 250    */
 251   z = y * L10EB;
 252   z += x * L10EB;
 253   z += e * L102B;
 254   z += y * L10EA;
 255   z += x * L10EA;
 256   z += e * L102A;
 257   return (z);
 258 }