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;
 186   double xhi;
 187
 188 /* Test for domain */
 189   xhi = ldbl_high (x);
 190   EXTRACT_WORDS64 (hx, xhi);
 191   if ((hx & 0x7fffffffffffffffLL) == 0)
 192     return (-1.0L / (x - x));
 193   if (hx < 0)
 194     return (x - x) / (x - x);
 195   if (hx >= 0x7ff0000000000000LL)
 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)