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