libquadmath/math/log1pq.c

   1 /*                                                      log1pl.c
   2  *
   3  *      Relative error logarithm
   4  *      Natural logarithm of 1+x for __float128 precision
   5  *
   6  *
   7  *
   8  * SYNOPSIS:
   9  *
  10  * __float128 x, y, log1pl();
  11  *
  12  * y = log1pq( x );
  13  *
  14  *
  15  *
  16  * DESCRIPTION:
  17  *
  18  * Returns the base e (2.718...) logarithm of 1+x.
  19  *
  20  * The argument 1+x is separated into its exponent and fractional
  21  * parts.  If the exponent is between -1 and +1, the logarithm
  22  * of the fraction is approximated by
  23  *
  24  *     log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
  25  *
  26  * Otherwise, setting  z = 2(w-1)/(w+1),
  27  *
  28  *     log(w) = z + z^3 P(z)/Q(z).
  29  *
  30  *
  31  *
  32  * ACCURACY:
  33  *
  34  *                      Relative error:
  35  * arithmetic   domain     # trials      peak         rms
  36  *    IEEE      -1, 8       100000      1.9e-34     4.3e-35
  37  */
  38
  39 /* Copyright 2001 by Stephen L. Moshier
  40
  41     This library is free software; you can redistribute it and/or
  42     modify it under the terms of the GNU Lesser General Public
  43     License as published by the Free Software Foundation; either
  44     version 2.1 of the License, or (at your option) any later version.
  45
  46     This library is distributed in the hope that it will be useful,
  47     but WITHOUT ANY WARRANTY; without even the implied warranty of
  48     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  49     Lesser General Public License for more details.
  50
  51     You should have received a copy of the GNU Lesser General Public
  52     License along with this library; if not, write to the Free Software
  53     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
  54
  55
  56 #include "quadmath-imp.h"
  57
  58 /* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
  59  * 1/sqrt(2) <= 1+x < sqrt(2)
  60  * Theoretical peak relative error = 5.3e-37,
  61  * relative peak error spread = 2.3e-14
  62  */
  63 static const __float128
  64   P12 = 1.538612243596254322971797716843006400388E-6Q,
  65   P11 = 4.998469661968096229986658302195402690910E-1Q,
  66   P10 = 2.321125933898420063925789532045674660756E1Q,
  67   P9 = 4.114517881637811823002128927449878962058E2Q,
  68   P8 = 3.824952356185897735160588078446136783779E3Q,
  69   P7 = 2.128857716871515081352991964243375186031E4Q,
  70   P6 = 7.594356839258970405033155585486712125861E4Q,
  71   P5 = 1.797628303815655343403735250238293741397E5Q,
  72   P4 = 2.854829159639697837788887080758954924001E5Q,
  73   P3 = 3.007007295140399532324943111654767187848E5Q,
  74   P2 = 2.014652742082537582487669938141683759923E5Q,
  75   P1 = 7.771154681358524243729929227226708890930E4Q,
  76   P0 = 1.313572404063446165910279910527789794488E4Q,
  77   /* Q12 = 1.000000000000000000000000000000000000000E0Q, */
  78   Q11 = 4.839208193348159620282142911143429644326E1Q,
  79   Q10 = 9.104928120962988414618126155557301584078E2Q,
  80   Q9 = 9.147150349299596453976674231612674085381E3Q,
  81   Q8 = 5.605842085972455027590989944010492125825E4Q,
  82   Q7 = 2.248234257620569139969141618556349415120E5Q,
  83   Q6 = 6.132189329546557743179177159925690841200E5Q,
  84   Q5 = 1.158019977462989115839826904108208787040E6Q,
  85   Q4 = 1.514882452993549494932585972882995548426E6Q,
  86   Q3 = 1.347518538384329112529391120390701166528E6Q,
  87   Q2 = 7.777690340007566932935753241556479363645E5Q,
  88   Q1 = 2.626900195321832660448791748036714883242E5Q,
  89   Q0 = 3.940717212190338497730839731583397586124E4Q;
  90
  91 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
  92  * where z = 2(x-1)/(x+1)
  93  * 1/sqrt(2) <= x < sqrt(2)
  94  * Theoretical peak relative error = 1.1e-35,
  95  * relative peak error spread 1.1e-9
  96  */
  97 static const __float128
  98   R5 = -8.828896441624934385266096344596648080902E-1Q,
  99   R4 = 8.057002716646055371965756206836056074715E1Q,
 100   R3 = -2.024301798136027039250415126250455056397E3Q,
 101   R2 = 2.048819892795278657810231591630928516206E4Q,
 102   R1 = -8.977257995689735303686582344659576526998E4Q,
 103   R0 = 1.418134209872192732479751274970992665513E5Q,
 104   /* S6 = 1.000000000000000000000000000000000000000E0Q, */
 105   S5 = -1.186359407982897997337150403816839480438E2Q,
 106   S4 = 3.998526750980007367835804959888064681098E3Q,
 107   S3 = -5.748542087379434595104154610899551484314E4Q,
 108   S2 = 4.001557694070773974936904547424676279307E5Q,
 109   S1 = -1.332535117259762928288745111081235577029E6Q,
 110   S0 = 1.701761051846631278975701529965589676574E6Q;
 111
 112 /* C1 + C2 = ln 2 */
 113 static const __float128 C1 = 6.93145751953125E-1Q;
 114 static const __float128 C2 = 1.428606820309417232121458176568075500134E-6Q;
 115
 116 static const __float128 sqrth = 0.7071067811865475244008443621048490392848Q;
 117 static const __float128 zero = 0.0Q;
 118
 119
 120 __float128
 121 log1pq (__float128 xm1)
 122 {
 123   __float128 x, y, z, r, s;
 124   ieee854_float128 u;
 125   int32_t hx;
 126   int e;
 127
 128   /* Test for NaN or infinity input. */
 129   u.value = xm1;
 130   hx = u.words32.w0;
 131   if (hx >= 0x7fff0000)
 132     return xm1;
 133
 134   /* log1p(+- 0) = +- 0.  */
 135   if (((hx & 0x7fffffff) == 0)
 136       && (u.words32.w1 | u.words32.w2 | u.words32.w3) == 0)
 137     return xm1;
 138
 139   if ((hx & 0x7fffffff) < 0x3f8e0000)
 140     {
 141       if ((int) xm1 == 0)
 142        return xm1;
 143     }
 144
 145   x = xm1 + 1.0Q;
 146
 147   /* log1p(-1) = -inf */
 148   if (x <= 0.0Q)
 149     {
 150       if (x == 0.0Q)
 151         return (-1.0Q / (x - x));
 152       else
 153         return (zero / (x - x));
 154     }
 155
 156   /* Separate mantissa from exponent.  */
 157
 158   /* Use frexp used so that denormal numbers will be handled properly.  */
 159   x = frexpq (x, &e);
 160
 161   /* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
 162      where z = 2(x-1)/x+1).  */
 163   if ((e > 2) || (e < -2))
 164     {
 165       if (x < sqrth)
 166         {                       /* 2( 2x-1 )/( 2x+1 ) */
 167           e -= 1;
 168           z = x - 0.5Q;
 169           y = 0.5Q * z + 0.5Q;
 170         }
 171       else
 172         {                       /*  2 (x-1)/(x+1)   */
 173           z = x - 0.5Q;
 174           z -= 0.5Q;
 175           y = 0.5Q * x + 0.5Q;
 176         }
 177       x = z / y;
 178       z = x * x;
 179       r = ((((R5 * z
 180               + R4) * z
 181              + R3) * z
 182             + R2) * z
 183            + R1) * z
 184         + R0;
 185       s = (((((z
 186                + S5) * z
 187               + S4) * z
 188              + S3) * z
 189             + S2) * z
 190            + S1) * z
 191         + S0;
 192       z = x * (z * r / s);
 193       z = z + e * C2;
 194       z = z + x;
 195       z = z + e * C1;
 196       return (z);
 197     }
 198
 199
 200   /* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
 201
 202   if (x < sqrth)
 203     {
 204       e -= 1;
 205       if (e != 0)
 206         x = 2.0Q * x - 1.0Q;    /*  2x - 1  */
 207       else
 208         x = xm1;
 209     }
 210   else
 211     {
 212       if (e != 0)
 213         x = x - 1.0Q;
 214       else
 215         x = xm1;
 216     }
 217   z = x * x;
 218   r = (((((((((((P12 * x
 219                  + P11) * x
 220                 + P10) * x
 221                + P9) * x
 222               + P8) * x
 223              + P7) * x
 224             + P6) * x
 225            + P5) * x
 226           + P4) * x
 227          + P3) * x
 228         + P2) * x
 229        + P1) * x
 230     + P0;
 231   s = (((((((((((x
 232                  + Q11) * x
 233                 + Q10) * x
 234                + Q9) * x
 235               + Q8) * x
 236              + Q7) * x
 237             + Q6) * x
 238            + Q5) * x
 239           + Q4) * x
 240          + Q3) * x
 241         + Q2) * x
 242        + Q1) * x
 243     + Q0;
 244   y = x * (z * r / s);
 245   y = y + e * C2;
 246   z = y - 0.5Q * z;
 247   z = z + x;
 248   z = z + e * C1;
 249   return (z);
 250 }