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 long double P[13] =
  71 {
  72   1.313572404063446165910279910527789794488E4L,
  73   7.771154681358524243729929227226708890930E4L,
  74   2.014652742082537582487669938141683759923E5L,
  75   3.007007295140399532324943111654767187848E5L,
  76   2.854829159639697837788887080758954924001E5L,
  77   1.797628303815655343403735250238293741397E5L,
  78   7.594356839258970405033155585486712125861E4L,
  79   2.128857716871515081352991964243375186031E4L,
  80   3.824952356185897735160588078446136783779E3L,
  81   4.114517881637811823002128927449878962058E2L,
  82   2.321125933898420063925789532045674660756E1L,
  83   4.998469661968096229986658302195402690910E-1L,
  84   1.538612243596254322971797716843006400388E-6L
  85 };
  86 static const long double Q[12] =
  87 {
  88   3.940717212190338497730839731583397586124E4L,
  89   2.626900195321832660448791748036714883242E5L,
  90   7.777690340007566932935753241556479363645E5L,
  91   1.347518538384329112529391120390701166528E6L,
  92   1.514882452993549494932585972882995548426E6L,
  93   1.158019977462989115839826904108208787040E6L,
  94   6.132189329546557743179177159925690841200E5L,
  95   2.248234257620569139969141618556349415120E5L,
  96   5.605842085972455027590989944010492125825E4L,
  97   9.147150349299596453976674231612674085381E3L,
  98   9.104928120962988414618126155557301584078E2L,
  99   4.839208193348159620282142911143429644326E1L
 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 long double R[6] =
 110 {
 111   1.418134209872192732479751274970992665513E5L,
 112  -8.977257995689735303686582344659576526998E4L,
 113   2.048819892795278657810231591630928516206E4L,
 114  -2.024301798136027039250415126250455056397E3L,
 115   8.057002716646055371965756206836056074715E1L,
 116  -8.828896441624934385266096344596648080902E-1L
 117 };
 118 static const long double S[6] =
 119 {
 120   1.701761051846631278975701529965589676574E6L,
 121  -1.332535117259762928288745111081235577029E6L,
 122   4.001557694070773974936904547424676279307E5L,
 123  -5.748542087379434595104154610899551484314E4L,
 124   3.998526750980007367835804959888064681098E3L,
 125  -1.186359407982897997337150403816839480438E2L
 126 /* 1.000000000000000000000000000000000000000E0L, */
 127 };
 128
 129 static const long double
 130 /* log2(e) - 1 */
 131 LOG2EA = 4.4269504088896340735992468100189213742664595E-1L,
 132 /* sqrt(2)/2 */
 133 SQRTH = 7.071067811865475244008443621048490392848359E-1L;
 134
 135
 136 /* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 137
 138 static long double
 139 neval (long double x, const long double *p, int n)
 140 {
 141   long double 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 long double
 157 deval (long double x, const long double *p, int n)
 158 {
 159   long double 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 long double
 174 __ieee754_log2l (x)
 175      long double x;
 176 {
 177   long double z;
 178   long double 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.0L / (x - x));
 186   if (hx < 0)
 187     return (x - x) / (x - x);
 188   if (hx >= 0x7fff000000000000LL)
 189     return (x + x);
 190
 191 /* separate mantissa from exponent */
 192
 193 /* Note, frexp is used so that denormal numbers
 194  * will be handled properly.
 195  */
 196   x = __frexpl (x, &e);
 197
 198
 199 /* logarithm using log(x) = z + z**3 P(z)/Q(z),
 200  * where z = 2(x-1)/x+1)
 201  */
 202   if ((e > 2) || (e < -2))
 203     {
 204       if (x < SQRTH)
 205         {                       /* 2( 2x-1 )/( 2x+1 ) */
 206           e -= 1;
 207           z = x - 0.5L;
 208           y = 0.5L * z + 0.5L;
 209         }
 210       else
 211         {                       /*  2 (x-1)/(x+1)   */
 212           z = x - 0.5L;
 213           z -= 0.5L;
 214           y = 0.5L * x + 0.5L;
 215         }
 216       x = z / y;
 217       z = x * x;
 218       y = x * (z * neval (z, R, 5) / deval (z, S, 5));
 219       goto done;
 220     }
 221
 222
 223 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
 224
 225   if (x < SQRTH)
 226     {
 227       e -= 1;
 228       x = 2.0 * x - 1.0L;       /*  2x - 1  */
 229     }
 230   else
 231     {
 232       x = x - 1.0L;
 233     }
 234   z = x * x;
 235   y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
 236   y = y - 0.5 * z;
 237
 238 done:
 239
 240 /* Multiply log of fraction by log2(e)
 241  * and base 2 exponent by 1
 242  */
 243   z = y * LOG2EA;
 244   z += x * LOG2EA;
 245   z += y;
 246   z += x;
 247   z += e;
 248   return (z);
 249 }
 250 strong_alias (__ieee754_log2l, __log2l_finite)