libquadmath/math/log10q.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 "quadmath-imp.h"
  64
  65 /* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
  66  * 1/sqrt(2) <= x < sqrt(2)
  67  * Theoretical peak relative error = 5.3e-37,
  68  * relative peak error spread = 2.3e-14
  69  */
  70 static const __float128 P[13] =
  71 {
  72   1.313572404063446165910279910527789794488E4Q,
  73   7.771154681358524243729929227226708890930E4Q,
  74   2.014652742082537582487669938141683759923E5Q,
  75   3.007007295140399532324943111654767187848E5Q,
  76   2.854829159639697837788887080758954924001E5Q,
  77   1.797628303815655343403735250238293741397E5Q,
  78   7.594356839258970405033155585486712125861E4Q,
  79   2.128857716871515081352991964243375186031E4Q,
  80   3.824952356185897735160588078446136783779E3Q,
  81   4.114517881637811823002128927449878962058E2Q,
  82   2.321125933898420063925789532045674660756E1Q,
  83   4.998469661968096229986658302195402690910E-1Q,
  84   1.538612243596254322971797716843006400388E-6Q
  85 };
  86 static const __float128 Q[12] =
  87 {
  88   3.940717212190338497730839731583397586124E4Q,
  89   2.626900195321832660448791748036714883242E5Q,
  90   7.777690340007566932935753241556479363645E5Q,
  91   1.347518538384329112529391120390701166528E6Q,
  92   1.514882452993549494932585972882995548426E6Q,
  93   1.158019977462989115839826904108208787040E6Q,
  94   6.132189329546557743179177159925690841200E5Q,
  95   2.248234257620569139969141618556349415120E5Q,
  96   5.605842085972455027590989944010492125825E4Q,
  97   9.147150349299596453976674231612674085381E3Q,
  98   9.104928120962988414618126155557301584078E2Q,
  99   4.839208193348159620282142911143429644326E1Q
 100 /* 1.000000000000000000000000000000000000000E0L, */
 101 };
 102
 103 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
 104  * where z = 2(x-1)/(x+1)
 105  * 1/sqrt(2) <= x < sqrt(2)
 106  * Theoretical peak relative error = 1.1e-35,
 107  * relative peak error spread 1.1e-9
 108  */
 109 static const __float128 R[6] =
 110 {
 111   1.418134209872192732479751274970992665513E5Q,
 112  -8.977257995689735303686582344659576526998E4Q,
 113   2.048819892795278657810231591630928516206E4Q,
 114  -2.024301798136027039250415126250455056397E3Q,
 115   8.057002716646055371965756206836056074715E1Q,
 116  -8.828896441624934385266096344596648080902E-1Q
 117 };
 118 static const __float128 S[6] =
 119 {
 120   1.701761051846631278975701529965589676574E6Q,
 121  -1.332535117259762928288745111081235577029E6Q,
 122   4.001557694070773974936904547424676279307E5Q,
 123  -5.748542087379434595104154610899551484314E4Q,
 124   3.998526750980007367835804959888064681098E3Q,
 125  -1.186359407982897997337150403816839480438E2Q
 126 /* 1.000000000000000000000000000000000000000E0L, */
 127 };
 128
 129 static const __float128
 130 /* log10(2) */
 131 L102A = 0.3125Q,
 132 L102B = -1.14700043360188047862611052755069732318101185E-2Q,
 133 /* log10(e) */
 134 L10EA = 0.5Q,
 135 L10EB = -6.570551809674817234887108108339491770560299E-2Q,
 136 /* sqrt(2)/2 */
 137 SQRTH = 7.071067811865475244008443621048490392848359E-1Q;
 138
 139
 140
 141 /* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 142
 143 static __float128
 144 neval (__float128 x, const __float128 *p, int n)
 145 {
 146   __float128 y;
 147
 148   p += n;
 149   y = *p--;
 150   do
 151     {
 152       y = y * x + *p--;
 153     }
 154   while (--n > 0);
 155   return y;
 156 }
 157
 158
 159 /* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 160
 161 static __float128
 162 deval (__float128 x, const __float128 *p, int n)
 163 {
 164   __float128 y;
 165
 166   p += n;
 167   y = x + *p--;
 168   do
 169     {
 170       y = y * x + *p--;
 171     }
 172   while (--n > 0);
 173   return y;
 174 }
 175
 176
 177
 178 __float128
 179 log10q (__float128 x)
 180 {
 181   __float128 z;
 182   __float128 y;
 183   int e;
 184   int64_t hx, lx;
 185
 186 /* Test for domain */
 187   GET_FLT128_WORDS64 (hx, lx, x);
 188   if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
 189     return (-1 / fabsq (x));            /* log10l(+-0)=-inf  */
 190   if (hx < 0)
 191     return (x - x) / (x - x);
 192   if (hx >= 0x7fff000000000000LL)
 193     return (x + x);
 194
 195   if (x == 1)
 196     return 0;
 197
 198 /* separate mantissa from exponent */
 199
 200 /* Note, frexp is used so that denormal numbers
 201  * will be handled properly.
 202  */
 203   x = frexpq (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.5Q;
 215           y = 0.5Q * z + 0.5Q;
 216         }
 217       else
 218         {                       /*  2 (x-1)/(x+1)   */
 219           z = x - 0.5Q;
 220           z -= 0.5Q;
 221           y = 0.5Q * x + 0.5Q;
 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;  /*  2x - 1  */
 236     }
 237   else
 238     {
 239       x = x - 1;
 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 }