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 <https://www.gnu.org/licenses/>.
  60  */
  61
  62 #include <math.h>
  63 #include <math_private.h>
  64 #include <libm-alias-finite.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 _Float128 P[13] =
  72 {
  73   L(1.313572404063446165910279910527789794488E4),
  74   L(7.771154681358524243729929227226708890930E4),
  75   L(2.014652742082537582487669938141683759923E5),
  76   L(3.007007295140399532324943111654767187848E5),
  77   L(2.854829159639697837788887080758954924001E5),
  78   L(1.797628303815655343403735250238293741397E5),
  79   L(7.594356839258970405033155585486712125861E4),
  80   L(2.128857716871515081352991964243375186031E4),
  81   L(3.824952356185897735160588078446136783779E3),
  82   L(4.114517881637811823002128927449878962058E2),
  83   L(2.321125933898420063925789532045674660756E1),
  84   L(4.998469661968096229986658302195402690910E-1),
  85   L(1.538612243596254322971797716843006400388E-6)
  86 };
  87 static const _Float128 Q[12] =
  88 {
  89   L(3.940717212190338497730839731583397586124E4),
  90   L(2.626900195321832660448791748036714883242E5),
  91   L(7.777690340007566932935753241556479363645E5),
  92   L(1.347518538384329112529391120390701166528E6),
  93   L(1.514882452993549494932585972882995548426E6),
  94   L(1.158019977462989115839826904108208787040E6),
  95   L(6.132189329546557743179177159925690841200E5),
  96   L(2.248234257620569139969141618556349415120E5),
  97   L(5.605842085972455027590989944010492125825E4),
  98   L(9.147150349299596453976674231612674085381E3),
  99   L(9.104928120962988414618126155557301584078E2),
 100   L(4.839208193348159620282142911143429644326E1)
 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 _Float128 R[6] =
 111 {
 112   L(1.418134209872192732479751274970992665513E5),
 113  L(-8.977257995689735303686582344659576526998E4),
 114   L(2.048819892795278657810231591630928516206E4),
 115  L(-2.024301798136027039250415126250455056397E3),
 116   L(8.057002716646055371965756206836056074715E1),
 117  L(-8.828896441624934385266096344596648080902E-1)
 118 };
 119 static const _Float128 S[6] =
 120 {
 121   L(1.701761051846631278975701529965589676574E6),
 122  L(-1.332535117259762928288745111081235577029E6),
 123   L(4.001557694070773974936904547424676279307E5),
 124  L(-5.748542087379434595104154610899551484314E4),
 125   L(3.998526750980007367835804959888064681098E3),
 126  L(-1.186359407982897997337150403816839480438E2)
 127 /* 1.000000000000000000000000000000000000000E0L, */
 128 };
 129
 130 static const _Float128
 131 /* log2(e) - 1 */
 132 LOG2EA = L(4.4269504088896340735992468100189213742664595E-1),
 133 /* sqrt(2)/2 */
 134 SQRTH = L(7.071067811865475244008443621048490392848359E-1);
 135
 136
 137 /* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 138
 139 static _Float128
 140 neval (_Float128 x, const _Float128 *p, int n)
 141 {
 142   _Float128 y;
 143
 144   p += n;
 145   y = *p--;
 146   do
 147     {
 148       y = y * x + *p--;
 149     }
 150   while (--n > 0);
 151   return y;
 152 }
 153
 154
 155 /* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 156
 157 static _Float128
 158 deval (_Float128 x, const _Float128 *p, int n)
 159 {
 160   _Float128 y;
 161
 162   p += n;
 163   y = x + *p--;
 164   do
 165     {
 166       y = y * x + *p--;
 167     }
 168   while (--n > 0);
 169   return y;
 170 }
 171
 172
 173
 174 _Float128
 175 __ieee754_log2l (_Float128 x)
 176 {
 177   _Float128 z;
 178   _Float128 y;
 179   int e;
 180   int64_t hx, lx;
 181
 182 /* Test for domain */
 183   GET_LDOUBLE_WORDS64 (hx, lx, x);
 184   if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
 185     return (-1 / fabsl (x));            /* log2l(+-0)=-inf  */
 186   if (hx < 0)
 187     return (x - x) / (x - x);
 188   if (hx >= 0x7fff000000000000LL)
 189     return (x + x);
 190
 191   if (x == 1)
 192     return 0;
 193
 194 /* separate mantissa from exponent */
 195
 196 /* Note, frexp is used so that denormal numbers
 197  * will be handled properly.
 198  */
 199   x = __frexpl (x, &e);
 200
 201
 202 /* logarithm using log(x) = z + z**3 P(z)/Q(z),
 203  * where z = 2(x-1)/x+1)
 204  */
 205   if ((e > 2) || (e < -2))
 206     {
 207       if (x < SQRTH)
 208         {                       /* 2( 2x-1 )/( 2x+1 ) */
 209           e -= 1;
 210           z = x - L(0.5);
 211           y = L(0.5) * z + L(0.5);
 212         }
 213       else
 214         {                       /*  2 (x-1)/(x+1)   */
 215           z = x - L(0.5);
 216           z -= L(0.5);
 217           y = L(0.5) * x + L(0.5);
 218         }
 219       x = z / y;
 220       z = x * x;
 221       y = x * (z * neval (z, R, 5) / deval (z, S, 5));
 222       goto done;
 223     }
 224
 225
 226 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
 227
 228   if (x < SQRTH)
 229     {
 230       e -= 1;
 231       x = 2.0 * x - 1;  /*  2x - 1  */
 232     }
 233   else
 234     {
 235       x = x - 1;
 236     }
 237   z = x * x;
 238   y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
 239   y = y - 0.5 * z;
 240
 241 done:
 242
 243 /* Multiply log of fraction by log2(e)
 244  * and base 2 exponent by 1
 245  */
 246   z = y * LOG2EA;
 247   z += x * LOG2EA;
 248   z += y;
 249   z += x;
 250   z += e;
 251   return (z);
 252 }
 253 libm_alias_finite (__ieee754_log2l, __log2l)