Back to index

glibc  2.9
s_log1pl.c
Go to the documentation of this file.
00001 /*                                               log1pl.c
00002  *
00003  *      Relative error logarithm
00004  *     Natural logarithm of 1+x, 128-bit long double precision
00005  *
00006  *
00007  *
00008  * SYNOPSIS:
00009  *
00010  * long double x, y, log1pl();
00011  *
00012  * y = log1pl( x );
00013  *
00014  *
00015  *
00016  * DESCRIPTION:
00017  *
00018  * Returns the base e (2.718...) logarithm of 1+x.
00019  *
00020  * The argument 1+x is separated into its exponent and fractional
00021  * parts.  If the exponent is between -1 and +1, the logarithm
00022  * of the fraction is approximated by
00023  *
00024  *     log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
00025  *
00026  * Otherwise, setting  z = 2(w-1)/(w+1),
00027  *
00028  *     log(w) = z + z^3 P(z)/Q(z).
00029  *
00030  *
00031  *
00032  * ACCURACY:
00033  *
00034  *                      Relative error:
00035  * arithmetic   domain     # trials      peak         rms
00036  *    IEEE      -1, 8       100000      1.9e-34     4.3e-35
00037  */
00038 
00039 /* Copyright 2001 by Stephen L. Moshier
00040 
00041     This library is free software; you can redistribute it and/or
00042     modify it under the terms of the GNU Lesser General Public
00043     License as published by the Free Software Foundation; either
00044     version 2.1 of the License, or (at your option) any later version.
00045 
00046     This library is distributed in the hope that it will be useful,
00047     but WITHOUT ANY WARRANTY; without even the implied warranty of
00048     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00049     Lesser General Public License for more details.
00050 
00051     You should have received a copy of the GNU Lesser General Public
00052     License along with this library; if not, write to the Free Software
00053     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
00054 
00055 
00056 #include "math.h"
00057 #include "math_private.h"
00058 #include <math_ldbl_opt.h>
00059 
00060 /* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
00061  * 1/sqrt(2) <= 1+x < sqrt(2)
00062  * Theoretical peak relative error = 5.3e-37,
00063  * relative peak error spread = 2.3e-14
00064  */
00065 static const long double
00066   P12 = 1.538612243596254322971797716843006400388E-6L,
00067   P11 = 4.998469661968096229986658302195402690910E-1L,
00068   P10 = 2.321125933898420063925789532045674660756E1L,
00069   P9 = 4.114517881637811823002128927449878962058E2L,
00070   P8 = 3.824952356185897735160588078446136783779E3L,
00071   P7 = 2.128857716871515081352991964243375186031E4L,
00072   P6 = 7.594356839258970405033155585486712125861E4L,
00073   P5 = 1.797628303815655343403735250238293741397E5L,
00074   P4 = 2.854829159639697837788887080758954924001E5L,
00075   P3 = 3.007007295140399532324943111654767187848E5L,
00076   P2 = 2.014652742082537582487669938141683759923E5L,
00077   P1 = 7.771154681358524243729929227226708890930E4L,
00078   P0 = 1.313572404063446165910279910527789794488E4L,
00079   /* Q12 = 1.000000000000000000000000000000000000000E0L, */
00080   Q11 = 4.839208193348159620282142911143429644326E1L,
00081   Q10 = 9.104928120962988414618126155557301584078E2L,
00082   Q9 = 9.147150349299596453976674231612674085381E3L,
00083   Q8 = 5.605842085972455027590989944010492125825E4L,
00084   Q7 = 2.248234257620569139969141618556349415120E5L,
00085   Q6 = 6.132189329546557743179177159925690841200E5L,
00086   Q5 = 1.158019977462989115839826904108208787040E6L,
00087   Q4 = 1.514882452993549494932585972882995548426E6L,
00088   Q3 = 1.347518538384329112529391120390701166528E6L,
00089   Q2 = 7.777690340007566932935753241556479363645E5L,
00090   Q1 = 2.626900195321832660448791748036714883242E5L,
00091   Q0 = 3.940717212190338497730839731583397586124E4L;
00092 
00093 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
00094  * where z = 2(x-1)/(x+1)
00095  * 1/sqrt(2) <= x < sqrt(2)
00096  * Theoretical peak relative error = 1.1e-35,
00097  * relative peak error spread 1.1e-9
00098  */
00099 static const long double
00100   R5 = -8.828896441624934385266096344596648080902E-1L,
00101   R4 = 8.057002716646055371965756206836056074715E1L,
00102   R3 = -2.024301798136027039250415126250455056397E3L,
00103   R2 = 2.048819892795278657810231591630928516206E4L,
00104   R1 = -8.977257995689735303686582344659576526998E4L,
00105   R0 = 1.418134209872192732479751274970992665513E5L,
00106   /* S6 = 1.000000000000000000000000000000000000000E0L, */
00107   S5 = -1.186359407982897997337150403816839480438E2L,
00108   S4 = 3.998526750980007367835804959888064681098E3L,
00109   S3 = -5.748542087379434595104154610899551484314E4L,
00110   S2 = 4.001557694070773974936904547424676279307E5L,
00111   S1 = -1.332535117259762928288745111081235577029E6L,
00112   S0 = 1.701761051846631278975701529965589676574E6L;
00113 
00114 /* C1 + C2 = ln 2 */
00115 static const long double C1 = 6.93145751953125E-1L;
00116 static const long double C2 = 1.428606820309417232121458176568075500134E-6L;
00117 
00118 static const long double sqrth = 0.7071067811865475244008443621048490392848L;
00119 /* ln (2^16384 * (1 - 2^-113)) */
00120 static const long double maxlog = 1.1356523406294143949491931077970764891253E4L;
00121 static const long double big = 2e300L;
00122 static const long double zero = 0.0L;
00123 
00124 
00125 long double
00126 __log1pl (long double xm1)
00127 {
00128   long double x, y, z, r, s;
00129   ieee854_long_double_shape_type u;
00130   int32_t hx;
00131   int e;
00132 
00133   /* Test for NaN or infinity input. */
00134   u.value = xm1;
00135   hx = u.parts32.w0;
00136   if (hx >= 0x7ff00000)
00137     return xm1;
00138 
00139   /* log1p(+- 0) = +- 0.  */
00140   if (((hx & 0x7fffffff) == 0)
00141       && (u.parts32.w1 | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3) == 0)
00142     return xm1;
00143 
00144   x = xm1 + 1.0L;
00145 
00146   /* log1p(-1) = -inf */
00147   if (x <= 0.0L)
00148     {
00149       if (x == 0.0L)
00150        return (-1.0L / (x - x));
00151       else
00152        return (zero / (x - x));
00153     }
00154 
00155   /* Separate mantissa from exponent.  */
00156 
00157   /* Use frexp used so that denormal numbers will be handled properly.  */
00158   x = __frexpl (x, &e);
00159 
00160   /* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
00161      where z = 2(x-1)/x+1).  */
00162   if ((e > 2) || (e < -2))
00163     {
00164       if (x < sqrth)
00165        {                    /* 2( 2x-1 )/( 2x+1 ) */
00166          e -= 1;
00167          z = x - 0.5L;
00168          y = 0.5L * z + 0.5L;
00169        }
00170       else
00171        {                    /*  2 (x-1)/(x+1)   */
00172          z = x - 0.5L;
00173          z -= 0.5L;
00174          y = 0.5L * x + 0.5L;
00175        }
00176       x = z / y;
00177       z = x * x;
00178       r = ((((R5 * z
00179              + R4) * z
00180             + R3) * z
00181            + R2) * z
00182           + R1) * z
00183        + R0;
00184       s = (((((z
00185               + S5) * z
00186              + S4) * z
00187             + S3) * z
00188            + S2) * z
00189           + S1) * z
00190        + S0;
00191       z = x * (z * r / s);
00192       z = z + e * C2;
00193       z = z + x;
00194       z = z + e * C1;
00195       return (z);
00196     }
00197 
00198 
00199   /* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
00200 
00201   if (x < sqrth)
00202     {
00203       e -= 1;
00204       if (e != 0)
00205        x = 2.0L * x - 1.0L; /*  2x - 1  */
00206       else
00207        x = xm1;
00208     }
00209   else
00210     {
00211       if (e != 0)
00212        x = x - 1.0L;
00213       else
00214        x = xm1;
00215     }
00216   z = x * x;
00217   r = (((((((((((P12 * x
00218                + P11) * x
00219               + P10) * x
00220               + P9) * x
00221              + P8) * x
00222             + P7) * x
00223            + P6) * x
00224           + P5) * x
00225          + P4) * x
00226         + P3) * x
00227        + P2) * x
00228        + P1) * x
00229     + P0;
00230   s = (((((((((((x
00231                + Q11) * x
00232               + Q10) * x
00233               + Q9) * x
00234              + Q8) * x
00235             + Q7) * x
00236            + Q6) * x
00237           + Q5) * x
00238          + Q4) * x
00239         + Q3) * x
00240        + Q2) * x
00241        + Q1) * x
00242     + Q0;
00243   y = x * (z * r / s);
00244   y = y + e * C2;
00245   z = y - 0.5L * z;
00246   z = z + x;
00247   z = z + e * C1;
00248   return (z);
00249 }
00250 
00251 long_double_symbol (libm, __log1pl, log1pl);