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, write to the Free Software
  61     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA
  62
  63  */
  64
  65 #include "quadmath-imp.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 __float128 P[13] =
  73 {
  74   1.313572404063446165910279910527789794488E4Q,
  75   7.771154681358524243729929227226708890930E4Q,
  76   2.014652742082537582487669938141683759923E5Q,
  77   3.007007295140399532324943111654767187848E5Q,
  78   2.854829159639697837788887080758954924001E5Q,
  79   1.797628303815655343403735250238293741397E5Q,
  80   7.594356839258970405033155585486712125861E4Q,
  81   2.128857716871515081352991964243375186031E4Q,
  82   3.824952356185897735160588078446136783779E3Q,
  83   4.114517881637811823002128927449878962058E2Q,
  84   2.321125933898420063925789532045674660756E1Q,
  85   4.998469661968096229986658302195402690910E-1Q,
  86   1.538612243596254322971797716843006400388E-6Q
  87 };
  88 static const __float128 Q[12] =
  89 {
  90   3.940717212190338497730839731583397586124E4Q,
  91   2.626900195321832660448791748036714883242E5Q,
  92   7.777690340007566932935753241556479363645E5Q,
  93   1.347518538384329112529391120390701166528E6Q,
  94   1.514882452993549494932585972882995548426E6Q,
  95   1.158019977462989115839826904108208787040E6Q,
  96   6.132189329546557743179177159925690841200E5Q,
  97   2.248234257620569139969141618556349415120E5Q,
  98   5.605842085972455027590989944010492125825E4Q,
  99   9.147150349299596453976674231612674085381E3Q,
 100   9.104928120962988414618126155557301584078E2Q,
 101   4.839208193348159620282142911143429644326E1Q
 102 /* 1.000000000000000000000000000000000000000E0Q, */
 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 __float128 R[6] =
 112 {
 113   1.418134209872192732479751274970992665513E5Q,
 114  -8.977257995689735303686582344659576526998E4Q,
 115   2.048819892795278657810231591630928516206E4Q,
 116  -2.024301798136027039250415126250455056397E3Q,
 117   8.057002716646055371965756206836056074715E1Q,
 118  -8.828896441624934385266096344596648080902E-1Q
 119 };
 120 static const __float128 S[6] =
 121 {
 122   1.701761051846631278975701529965589676574E6Q,
 123  -1.332535117259762928288745111081235577029E6Q,
 124   4.001557694070773974936904547424676279307E5Q,
 125  -5.748542087379434595104154610899551484314E4Q,
 126   3.998526750980007367835804959888064681098E3Q,
 127  -1.186359407982897997337150403816839480438E2Q
 128 /* 1.000000000000000000000000000000000000000E0Q, */
 129 };
 130
 131 static const __float128
 132 /* log10(2) */
 133 L102A = 0.3125Q,
 134 L102B = -1.14700043360188047862611052755069732318101185E-2Q,
 135 /* log10(e) */
 136 L10EA = 0.5Q,
 137 L10EB = -6.570551809674817234887108108339491770560299E-2Q,
 138 /* sqrt(2)/2 */
 139 SQRTH = 7.071067811865475244008443621048490392848359E-1Q;
 140
 141
 142
 143 /* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 144
 145 static __float128
 146 neval (__float128 x, const __float128 *p, int n)
 147 {
 148   __float128 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 __float128
 164 deval (__float128 x, const __float128 *p, int n)
 165 {
 166   __float128 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 __float128
 181 log10q (__float128 x)
 182 {
 183   __float128 z;
 184   __float128 y;
 185   int e;
 186   int64_t hx, lx;
 187
 188 /* Test for domain */
 189   GET_FLT128_WORDS64 (hx, lx, x);
 190   if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
 191     return (-1.0Q / (x - x));
 192   if (hx < 0)
 193     return (x - x) / (x - x);
 194   if (hx >= 0x7fff000000000000LL)
 195     return (x + x);
 196
 197 /* separate mantissa from exponent */
 198
 199 /* Note, frexp is used so that denormal numbers
 200  * will be handled properly.
 201  */
 202   x = frexpq (x, &e);
 203
 204
 205 /* logarithm using log(x) = z + z**3 P(z)/Q(z),
 206  * where z = 2(x-1)/x+1)
 207  */
 208   if ((e > 2) || (e < -2))
 209     {
 210       if (x < SQRTH)
 211         {                       /* 2( 2x-1 )/( 2x+1 ) */
 212           e -= 1;
 213           z = x - 0.5Q;
 214           y = 0.5Q * z + 0.5Q;
 215         }
 216       else
 217         {                       /*  2 (x-1)/(x+1)   */
 218           z = x - 0.5Q;
 219           z -= 0.5Q;
 220           y = 0.5Q * x + 0.5Q;
 221         }
 222       x = z / y;
 223       z = x * x;
 224       y = x * (z * neval (z, R, 5) / deval (z, S, 5));
 225       goto done;
 226     }
 227
 228
 229 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
 230
 231   if (x < SQRTH)
 232     {
 233       e -= 1;
 234       x = 2.0 * x - 1.0Q;       /*  2x - 1  */
 235     }
 236   else
 237     {
 238       x = x - 1.0Q;
 239     }
 240   z = x * x;
 241   y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
 242   y = y - 0.5 * z;
 243
 244 done:
 245
 246   /* Multiply log of fraction by log10(e)
 247    * and base 2 exponent by log10(2).
 248    */
 249   z = y * L10EB;
 250   z += x * L10EB;
 251   z += e * L102B;
 252   z += y * L10EA;
 253   z += x * L10EA;
 254   z += e * L102A;
 255   return (z);
 256 }