libquadmath/math/log1pq.c

   1 /*                                                      log1pl.c
   2  *
   3  *      Relative error logarithm
   4  *      Natural logarithm of 1+x, 128-bit long double precision
   5  *
   6  *
   7  *
   8  * SYNOPSIS:
   9  *
  10  * long double x, y, log1pl();
  11  *
  12  * y = log1pl( 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   x = xm1 + 1.0Q;
 140
 141   /* log1p(-1) = -inf */
 142   if (x <= 0.0Q)
 143     {
 144       if (x == 0.0Q)
 145         return (-1.0Q / (x - x));
 146       else
 147         return (zero / (x - x));
 148     }
 149
 150   /* Separate mantissa from exponent.  */
 151
 152   /* Use frexp used so that denormal numbers will be handled properly.  */
 153   x = frexpq (x, &e);
 154
 155   /* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
 156      where z = 2(x-1)/x+1).  */
 157   if ((e > 2) || (e < -2))
 158     {
 159       if (x < sqrth)
 160         {                       /* 2( 2x-1 )/( 2x+1 ) */
 161           e -= 1;
 162           z = x - 0.5Q;
 163           y = 0.5Q * z + 0.5Q;
 164         }
 165       else
 166         {                       /*  2 (x-1)/(x+1)   */
 167           z = x - 0.5Q;
 168           z -= 0.5Q;
 169           y = 0.5Q * x + 0.5Q;
 170         }
 171       x = z / y;
 172       z = x * x;
 173       r = ((((R5 * z
 174               + R4) * z
 175              + R3) * z
 176             + R2) * z
 177            + R1) * z
 178         + R0;
 179       s = (((((z
 180                + S5) * z
 181               + S4) * z
 182              + S3) * z
 183             + S2) * z
 184            + S1) * z
 185         + S0;
 186       z = x * (z * r / s);
 187       z = z + e * C2;
 188       z = z + x;
 189       z = z + e * C1;
 190       return (z);
 191     }
 192
 193
 194   /* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
 195
 196   if (x < sqrth)
 197     {
 198       e -= 1;
 199       if (e != 0)
 200         x = 2.0Q * x - 1.0Q;    /*  2x - 1  */
 201       else
 202         x = xm1;
 203     }
 204   else
 205     {
 206       if (e != 0)
 207         x = x - 1.0Q;
 208       else
 209         x = xm1;
 210     }
 211   z = x * x;
 212   r = (((((((((((P12 * x
 213                  + P11) * x
 214                 + P10) * x
 215                + P9) * x
 216               + P8) * x
 217              + P7) * x
 218             + P6) * x
 219            + P5) * x
 220           + P4) * x
 221          + P3) * x
 222         + P2) * x
 223        + P1) * x
 224     + P0;
 225   s = (((((((((((x
 226                  + Q11) * x
 227                 + Q10) * x
 228                + Q9) * x
 229               + Q8) * x
 230              + Q7) * x
 231             + Q6) * x
 232            + Q5) * x
 233           + Q4) * x
 234          + Q3) * x
 235         + Q2) * x
 236        + Q1) * x
 237     + Q0;
 238   y = x * (z * r / s);
 239   y = y + e * C2;
 240   z = y - 0.5Q * z;
 241   z = z + x;
 242   z = z + e * C1;
 243   return (z);
 244 }