sysdeps/ieee754/ldbl-128/e_log2l.c

   1 /*                                                      log2l.c
   2  *      Base 2 logarithm, 128-bit long double precision
   3  *
   4  *
   5  *
   6  * SYNOPSIS:
   7  *
   8  * long double x, y, log2l();
   9  *
  10  * y = log2l( x );
  11  *
  12  *
  13  *
  14  * DESCRIPTION:
  15  *
  16  * Returns the base 2 logarithm of x.
  17  *
  18  * The argument is separated into its exponent and fractional
  19  * parts.  If the exponent is between -1 and +1, the (natural)
  20  * logarithm of the fraction is approximated by
  21  *
  22  *     log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
  23  *
  24  * Otherwise, setting  z = 2(x-1)/x+1),
  25  *
  26  *     log(x) = z + z^3 P(z)/Q(z).
  27  *
  28  *
  29  *
  30  * ACCURACY:
  31  *
  32  *                      Relative error:
  33  * arithmetic   domain     # trials      peak         rms
  34  *    IEEE      0.5, 2.0     100,000    2.6e-34     4.9e-35
  35  *    IEEE     exp(+-10000)  100,000    9.6e-35     4.0e-35
  36  *
  37  * In the tests over the interval exp(+-10000), the logarithms
  38  * of the random arguments were uniformly distributed over
  39  * [-10000, +10000].
  40  *
  41  */
  42
  43 /*
  44    Cephes Math Library Release 2.2:  January, 1991
  45    Copyright 1984, 1991 by Stephen L. Moshier
  46    Adapted for glibc November, 2001
  47
  48     This library is free software; you can redistribute it and/or
  49     modify it under the terms of the GNU Lesser General Public
  50     License as published by the Free Software Foundation; either
  51     version 2.1 of the License, or (at your option) any later version.
  52
  53     This library is distributed in the hope that it will be useful,
  54     but WITHOUT ANY WARRANTY; without even the implied warranty of
  55     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  56     Lesser General Public License for more details.
  57
  58     You should have received a copy of the GNU Lesser General Public
  59     License along with this library; if not, see <http://www.gnu.org/licenses/>.
  60  */
  61
  62 #include <math.h>
  63 #include <math_private.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   L(1.313572404063446165910279910527789794488E4),
  73   L(7.771154681358524243729929227226708890930E4),
  74   L(2.014652742082537582487669938141683759923E5),
  75   L(3.007007295140399532324943111654767187848E5),
  76   L(2.854829159639697837788887080758954924001E5),
  77   L(1.797628303815655343403735250238293741397E5),
  78   L(7.594356839258970405033155585486712125861E4),
  79   L(2.128857716871515081352991964243375186031E4),
  80   L(3.824952356185897735160588078446136783779E3),
  81   L(4.114517881637811823002128927449878962058E2),
  82   L(2.321125933898420063925789532045674660756E1),
  83   L(4.998469661968096229986658302195402690910E-1),
  84   L(1.538612243596254322971797716843006400388E-6)
  85 };
  86 static const _Float128 Q[12] =
  87 {
  88   L(3.940717212190338497730839731583397586124E4),
  89   L(2.626900195321832660448791748036714883242E5),
  90   L(7.777690340007566932935753241556479363645E5),
  91   L(1.347518538384329112529391120390701166528E6),
  92   L(1.514882452993549494932585972882995548426E6),
  93   L(1.158019977462989115839826904108208787040E6),
  94   L(6.132189329546557743179177159925690841200E5),
  95   L(2.248234257620569139969141618556349415120E5),
  96   L(5.605842085972455027590989944010492125825E4),
  97   L(9.147150349299596453976674231612674085381E3),
  98   L(9.104928120962988414618126155557301584078E2),
  99   L(4.839208193348159620282142911143429644326E1)
 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   L(1.418134209872192732479751274970992665513E5),
 112  L(-8.977257995689735303686582344659576526998E4),
 113   L(2.048819892795278657810231591630928516206E4),
 114  L(-2.024301798136027039250415126250455056397E3),
 115   L(8.057002716646055371965756206836056074715E1),
 116  L(-8.828896441624934385266096344596648080902E-1)
 117 };
 118 static const _Float128 S[6] =
 119 {
 120   L(1.701761051846631278975701529965589676574E6),
 121  L(-1.332535117259762928288745111081235577029E6),
 122   L(4.001557694070773974936904547424676279307E5),
 123  L(-5.748542087379434595104154610899551484314E4),
 124   L(3.998526750980007367835804959888064681098E3),
 125  L(-1.186359407982897997337150403816839480438E2)
 126 /* 1.000000000000000000000000000000000000000E0L, */
 127 };
 128
 129 static const _Float128
 130 /* log2(e) - 1 */
 131 LOG2EA = L(4.4269504088896340735992468100189213742664595E-1),
 132 /* sqrt(2)/2 */
 133 SQRTH = L(7.071067811865475244008443621048490392848359E-1);
 134
 135
 136 /* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 137
 138 static _Float128
 139 neval (_Float128 x, const _Float128 *p, int n)
 140 {
 141   _Float128 y;
 142
 143   p += n;
 144   y = *p--;
 145   do
 146     {
 147       y = y * x + *p--;
 148     }
 149   while (--n > 0);
 150   return y;
 151 }
 152
 153
 154 /* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 155
 156 static _Float128
 157 deval (_Float128 x, const _Float128 *p, int n)
 158 {
 159   _Float128 y;
 160
 161   p += n;
 162   y = x + *p--;
 163   do
 164     {
 165       y = y * x + *p--;
 166     }
 167   while (--n > 0);
 168   return y;
 169 }
 170
 171
 172
 173 _Float128
 174 __ieee754_log2l (_Float128 x)
 175 {
 176   _Float128 z;
 177   _Float128 y;
 178   int e;
 179   int64_t hx, lx;
 180
 181 /* Test for domain */
 182   GET_LDOUBLE_WORDS64 (hx, lx, x);
 183   if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
 184     return (-1 / fabsl (x));            /* log2l(+-0)=-inf  */
 185   if (hx < 0)
 186     return (x - x) / (x - x);
 187   if (hx >= 0x7fff000000000000LL)
 188     return (x + x);
 189
 190   if (x == 1)
 191     return 0;
 192
 193 /* separate mantissa from exponent */
 194
 195 /* Note, frexp is used so that denormal numbers
 196  * will be handled properly.
 197  */
 198   x = __frexpl (x, &e);
 199
 200
 201 /* logarithm using log(x) = z + z**3 P(z)/Q(z),
 202  * where z = 2(x-1)/x+1)
 203  */
 204   if ((e > 2) || (e < -2))
 205     {
 206       if (x < SQRTH)
 207         {                       /* 2( 2x-1 )/( 2x+1 ) */
 208           e -= 1;
 209           z = x - L(0.5);
 210           y = L(0.5) * z + L(0.5);
 211         }
 212       else
 213         {                       /*  2 (x-1)/(x+1)   */
 214           z = x - L(0.5);
 215           z -= L(0.5);
 216           y = L(0.5) * x + L(0.5);
 217         }
 218       x = z / y;
 219       z = x * x;
 220       y = x * (z * neval (z, R, 5) / deval (z, S, 5));
 221       goto done;
 222     }
 223
 224
 225 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
 226
 227   if (x < SQRTH)
 228     {
 229       e -= 1;
 230       x = 2.0 * x - 1;  /*  2x - 1  */
 231     }
 232   else
 233     {
 234       x = x - 1;
 235     }
 236   z = x * x;
 237   y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
 238   y = y - 0.5 * z;
 239
 240 done:
 241
 242 /* Multiply log of fraction by log2(e)
 243  * and base 2 exponent by 1
 244  */
 245   z = y * LOG2EA;
 246   z += x * LOG2EA;
 247   z += y;
 248   z += x;
 249   z += e;
 250   return (z);
 251 }
 252 strong_alias (__ieee754_log2l, __log2l_finite)