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