libquadmath/math/log2q.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 "quadmath-imp.h"
  63
  64 /* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
  65  * 1/sqrt(2) <= x < sqrt(2)
  66  * Theoretical peak relative error = 5.3e-37,
  67  * relative peak error spread = 2.3e-14
  68  */
  69 static const __float128 P[13] =
  70 {
  71   1.313572404063446165910279910527789794488E4Q,
  72   7.771154681358524243729929227226708890930E4Q,
  73   2.014652742082537582487669938141683759923E5Q,
  74   3.007007295140399532324943111654767187848E5Q,
  75   2.854829159639697837788887080758954924001E5Q,
  76   1.797628303815655343403735250238293741397E5Q,
  77   7.594356839258970405033155585486712125861E4Q,
  78   2.128857716871515081352991964243375186031E4Q,
  79   3.824952356185897735160588078446136783779E3Q,
  80   4.114517881637811823002128927449878962058E2Q,
  81   2.321125933898420063925789532045674660756E1Q,
  82   4.998469661968096229986658302195402690910E-1Q,
  83   1.538612243596254322971797716843006400388E-6Q
  84 };
  85 static const __float128 Q[12] =
  86 {
  87   3.940717212190338497730839731583397586124E4Q,
  88   2.626900195321832660448791748036714883242E5Q,
  89   7.777690340007566932935753241556479363645E5Q,
  90   1.347518538384329112529391120390701166528E6Q,
  91   1.514882452993549494932585972882995548426E6Q,
  92   1.158019977462989115839826904108208787040E6Q,
  93   6.132189329546557743179177159925690841200E5Q,
  94   2.248234257620569139969141618556349415120E5Q,
  95   5.605842085972455027590989944010492125825E4Q,
  96   9.147150349299596453976674231612674085381E3Q,
  97   9.104928120962988414618126155557301584078E2Q,
  98   4.839208193348159620282142911143429644326E1Q
  99 /* 1.000000000000000000000000000000000000000E0L, */
 100 };
 101
 102 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
 103  * where z = 2(x-1)/(x+1)
 104  * 1/sqrt(2) <= x < sqrt(2)
 105  * Theoretical peak relative error = 1.1e-35,
 106  * relative peak error spread 1.1e-9
 107  */
 108 static const __float128 R[6] =
 109 {
 110   1.418134209872192732479751274970992665513E5Q,
 111  -8.977257995689735303686582344659576526998E4Q,
 112   2.048819892795278657810231591630928516206E4Q,
 113  -2.024301798136027039250415126250455056397E3Q,
 114   8.057002716646055371965756206836056074715E1Q,
 115  -8.828896441624934385266096344596648080902E-1Q
 116 };
 117 static const __float128 S[6] =
 118 {
 119   1.701761051846631278975701529965589676574E6Q,
 120  -1.332535117259762928288745111081235577029E6Q,
 121   4.001557694070773974936904547424676279307E5Q,
 122  -5.748542087379434595104154610899551484314E4Q,
 123   3.998526750980007367835804959888064681098E3Q,
 124  -1.186359407982897997337150403816839480438E2Q
 125 /* 1.000000000000000000000000000000000000000E0L, */
 126 };
 127
 128 static const __float128
 129 /* log2(e) - 1 */
 130 LOG2EA = 4.4269504088896340735992468100189213742664595E-1Q,
 131 /* sqrt(2)/2 */
 132 SQRTH = 7.071067811865475244008443621048490392848359E-1Q;
 133
 134
 135 /* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 136
 137 static __float128
 138 neval (__float128 x, const __float128 *p, int n)
 139 {
 140   __float128 y;
 141
 142   p += n;
 143   y = *p--;
 144   do
 145     {
 146       y = y * x + *p--;
 147     }
 148   while (--n > 0);
 149   return y;
 150 }
 151
 152
 153 /* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 154
 155 static __float128
 156 deval (__float128 x, const __float128 *p, int n)
 157 {
 158   __float128 y;
 159
 160   p += n;
 161   y = x + *p--;
 162   do
 163     {
 164       y = y * x + *p--;
 165     }
 166   while (--n > 0);
 167   return y;
 168 }
 169
 170
 171
 172 __float128
 173 log2q (__float128 x)
 174 {
 175   __float128 z;
 176   __float128 y;
 177   int e;
 178   int64_t hx, lx;
 179
 180 /* Test for domain */
 181   GET_FLT128_WORDS64 (hx, lx, x);
 182   if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
 183     return (-1 / fabsq (x));            /* log2l(+-0)=-inf  */
 184   if (hx < 0)
 185     return (x - x) / (x - x);
 186   if (hx >= 0x7fff000000000000LL)
 187     return (x + x);
 188
 189   if (x == 1)
 190     return 0;
 191
 192 /* separate mantissa from exponent */
 193
 194 /* Note, frexp is used so that denormal numbers
 195  * will be handled properly.
 196  */
 197   x = frexpq (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.5Q;
 209           y = 0.5Q * z + 0.5Q;
 210         }
 211       else
 212         {                       /*  2 (x-1)/(x+1)   */
 213           z = x - 0.5Q;
 214           z -= 0.5Q;
 215           y = 0.5Q * x + 0.5Q;
 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;  /*  2x - 1  */
 230     }
 231   else
 232     {
 233       x = x - 1;
 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 }