libquadmath/math/log2q.c

   1 /*                                                      log2q.c
   2  *      Base 2 logarithm for __float128 precision
   3  *
   4  *
   5  *
   6  * SYNOPSIS:
   7  *
   8  * __float128 x, y, log2q();
   9  *
  10  * y = log2q( 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 "quadmath-imp.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   1.313572404063446165910279910527789794488E4Q,
  73   7.771154681358524243729929227226708890930E4Q,
  74   2.014652742082537582487669938141683759923E5Q,
  75   3.007007295140399532324943111654767187848E5Q,
  76   2.854829159639697837788887080758954924001E5Q,
  77   1.797628303815655343403735250238293741397E5Q,
  78   7.594356839258970405033155585486712125861E4Q,
  79   2.128857716871515081352991964243375186031E4Q,
  80   3.824952356185897735160588078446136783779E3Q,
  81   4.114517881637811823002128927449878962058E2Q,
  82   2.321125933898420063925789532045674660756E1Q,
  83   4.998469661968096229986658302195402690910E-1Q,
  84   1.538612243596254322971797716843006400388E-6Q
  85 };
  86 static const __float128 Q[12] =
  87 {
  88   3.940717212190338497730839731583397586124E4Q,
  89   2.626900195321832660448791748036714883242E5Q,
  90   7.777690340007566932935753241556479363645E5Q,
  91   1.347518538384329112529391120390701166528E6Q,
  92   1.514882452993549494932585972882995548426E6Q,
  93   1.158019977462989115839826904108208787040E6Q,
  94   6.132189329546557743179177159925690841200E5Q,
  95   2.248234257620569139969141618556349415120E5Q,
  96   5.605842085972455027590989944010492125825E4Q,
  97   9.147150349299596453976674231612674085381E3Q,
  98   9.104928120962988414618126155557301584078E2Q,
  99   4.839208193348159620282142911143429644326E1Q
 100 /* 1.000000000000000000000000000000000000000E0Q, */
 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   1.418134209872192732479751274970992665513E5Q,
 112  -8.977257995689735303686582344659576526998E4Q,
 113   2.048819892795278657810231591630928516206E4Q,
 114  -2.024301798136027039250415126250455056397E3Q,
 115   8.057002716646055371965756206836056074715E1Q,
 116  -8.828896441624934385266096344596648080902E-1Q
 117 };
 118 static const __float128 S[6] =
 119 {
 120   1.701761051846631278975701529965589676574E6Q,
 121  -1.332535117259762928288745111081235577029E6Q,
 122   4.001557694070773974936904547424676279307E5Q,
 123  -5.748542087379434595104154610899551484314E4Q,
 124   3.998526750980007367835804959888064681098E3Q,
 125  -1.186359407982897997337150403816839480438E2Q
 126 /* 1.000000000000000000000000000000000000000E0Q, */
 127 };
 128
 129 static const __float128
 130 /* log2(e) - 1 */
 131 LOG2EA = 4.4269504088896340735992468100189213742664595E-1Q,
 132 /* sqrt(2)/2 */
 133 SQRTH = 7.071067811865475244008443621048490392848359E-1Q;
 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 log2q (__float128 x)
 175 {
 176   __float128 z;
 177   __float128 y;
 178   int e;
 179   int64_t hx, lx;
 180
 181 /* Test for domain */
 182   GET_FLT128_WORDS64 (hx, lx, x);
 183   if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
 184     return (-1.0Q / fabsq (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.0Q)
 191     return 0.0Q;
 192
 193 /* separate mantissa from exponent */
 194
 195 /* Note, frexp is used so that denormal numbers
 196  * will be handled properly.
 197  */
 198   x = frexpq (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 - 0.5Q;
 210           y = 0.5Q * z + 0.5Q;
 211         }
 212       else
 213         {                       /*  2 (x-1)/(x+1)   */
 214           z = x - 0.5Q;
 215           z -= 0.5Q;
 216           y = 0.5Q * x + 0.5Q;
 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.0Q;       /*  2x - 1  */
 231     }
 232   else
 233     {
 234       x = x - 1.0Q;
 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 }