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 
00059 /* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
00060  * 1/sqrt(2) <= 1+x < sqrt(2)
00061  * Theoretical peak relative error = 5.3e-37,
00062  * relative peak error spread = 2.3e-14
00063  */
00064 static const long double
00065   P12 = 1.538612243596254322971797716843006400388E-6L,
00066   P11 = 4.998469661968096229986658302195402690910E-1L,
00067   P10 = 2.321125933898420063925789532045674660756E1L,
00068   P9 = 4.114517881637811823002128927449878962058E2L,
00069   P8 = 3.824952356185897735160588078446136783779E3L,
00070   P7 = 2.128857716871515081352991964243375186031E4L,
00071   P6 = 7.594356839258970405033155585486712125861E4L,
00072   P5 = 1.797628303815655343403735250238293741397E5L,
00073   P4 = 2.854829159639697837788887080758954924001E5L,
00074   P3 = 3.007007295140399532324943111654767187848E5L,
00075   P2 = 2.014652742082537582487669938141683759923E5L,
00076   P1 = 7.771154681358524243729929227226708890930E4L,
00077   P0 = 1.313572404063446165910279910527789794488E4L,
00078   /* Q12 = 1.000000000000000000000000000000000000000E0L, */
00079   Q11 = 4.839208193348159620282142911143429644326E1L,
00080   Q10 = 9.104928120962988414618126155557301584078E2L,
00081   Q9 = 9.147150349299596453976674231612674085381E3L,
00082   Q8 = 5.605842085972455027590989944010492125825E4L,
00083   Q7 = 2.248234257620569139969141618556349415120E5L,
00084   Q6 = 6.132189329546557743179177159925690841200E5L,
00085   Q5 = 1.158019977462989115839826904108208787040E6L,
00086   Q4 = 1.514882452993549494932585972882995548426E6L,
00087   Q3 = 1.347518538384329112529391120390701166528E6L,
00088   Q2 = 7.777690340007566932935753241556479363645E5L,
00089   Q1 = 2.626900195321832660448791748036714883242E5L,
00090   Q0 = 3.940717212190338497730839731583397586124E4L;
00091 
00092 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
00093  * where z = 2(x-1)/(x+1)
00094  * 1/sqrt(2) <= x < sqrt(2)
00095  * Theoretical peak relative error = 1.1e-35,
00096  * relative peak error spread 1.1e-9
00097  */
00098 static const long double
00099   R5 = -8.828896441624934385266096344596648080902E-1L,
00100   R4 = 8.057002716646055371965756206836056074715E1L,
00101   R3 = -2.024301798136027039250415126250455056397E3L,
00102   R2 = 2.048819892795278657810231591630928516206E4L,
00103   R1 = -8.977257995689735303686582344659576526998E4L,
00104   R0 = 1.418134209872192732479751274970992665513E5L,
00105   /* S6 = 1.000000000000000000000000000000000000000E0L, */
00106   S5 = -1.186359407982897997337150403816839480438E2L,
00107   S4 = 3.998526750980007367835804959888064681098E3L,
00108   S3 = -5.748542087379434595104154610899551484314E4L,
00109   S2 = 4.001557694070773974936904547424676279307E5L,
00110   S1 = -1.332535117259762928288745111081235577029E6L,
00111   S0 = 1.701761051846631278975701529965589676574E6L;
00112 
00113 /* C1 + C2 = ln 2 */
00114 static const long double C1 = 6.93145751953125E-1L;
00115 static const long double C2 = 1.428606820309417232121458176568075500134E-6L;
00116 
00117 static const long double sqrth = 0.7071067811865475244008443621048490392848L;
00118 /* ln (2^16384 * (1 - 2^-113)) */
00119 static const long double maxlog = 1.1356523406294143949491931077970764891253E4L;
00120 static const long double big = 2e4932L;
00121 static const long double zero = 0.0L;
00122 
00123 long double
00124 __log1pl (long double xm1)
00125 {
00126   long double x, y, z, r, s;
00127   ieee854_long_double_shape_type u;
00128   int32_t hx;
00129   int e;
00130 
00131   /* Test for NaN or infinity input. */
00132   u.value = xm1;
00133   hx = u.parts32.w0;
00134   if (hx >= 0x7fff0000)
00135     return xm1;
00136 
00137   /* log1p(+- 0) = +- 0.  */
00138   if (((hx & 0x7fffffff) == 0)
00139       && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
00140     return xm1;
00141 
00142   x = xm1 + 1.0L;
00143 
00144   /* log1p(-1) = -inf */
00145   if (x <= 0.0L)
00146     {
00147       if (x == 0.0L)
00148        return (-1.0L / (x - x));
00149       else
00150        return (zero / (x - x));
00151     }
00152 
00153   /* Separate mantissa from exponent.  */
00154 
00155   /* Use frexp used so that denormal numbers will be handled properly.  */
00156   x = __frexpl (x, &e);
00157 
00158   /* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
00159      where z = 2(x-1)/x+1).  */
00160   if ((e > 2) || (e < -2))
00161     {
00162       if (x < sqrth)
00163        {                    /* 2( 2x-1 )/( 2x+1 ) */
00164          e -= 1;
00165          z = x - 0.5L;
00166          y = 0.5L * z + 0.5L;
00167        }
00168       else
00169        {                    /*  2 (x-1)/(x+1)   */
00170          z = x - 0.5L;
00171          z -= 0.5L;
00172          y = 0.5L * x + 0.5L;
00173        }
00174       x = z / y;
00175       z = x * x;
00176       r = ((((R5 * z
00177              + R4) * z
00178             + R3) * z
00179            + R2) * z
00180           + R1) * z
00181        + R0;
00182       s = (((((z
00183               + S5) * z
00184              + S4) * z
00185             + S3) * z
00186            + S2) * z
00187           + S1) * z
00188        + S0;
00189       z = x * (z * r / s);
00190       z = z + e * C2;
00191       z = z + x;
00192       z = z + e * C1;
00193       return (z);
00194     }
00195 
00196 
00197   /* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
00198 
00199   if (x < sqrth)
00200     {
00201       e -= 1;
00202       if (e != 0)
00203        x = 2.0L * x - 1.0L; /*  2x - 1  */
00204       else
00205        x = xm1;
00206     }
00207   else
00208     {
00209       if (e != 0)
00210        x = x - 1.0L;
00211       else
00212        x = xm1;
00213     }
00214   z = x * x;
00215   r = (((((((((((P12 * x
00216                + P11) * x
00217               + P10) * x
00218               + P9) * x
00219              + P8) * x
00220             + P7) * x
00221            + P6) * x
00222           + P5) * x
00223          + P4) * x
00224         + P3) * x
00225        + P2) * x
00226        + P1) * x
00227     + P0;
00228   s = (((((((((((x
00229                + Q11) * x
00230               + Q10) * x
00231               + Q9) * x
00232              + Q8) * x
00233             + Q7) * x
00234            + Q6) * x
00235           + Q5) * x
00236          + Q4) * x
00237         + Q3) * x
00238        + Q2) * x
00239        + Q1) * x
00240     + Q0;
00241   y = x * (z * r / s);
00242   y = y + e * C2;
00243   z = y - 0.5L * z;
00244   z = z + x;
00245   z = z + e * C1;
00246   return (z);
00247 }
00248 
00249 weak_alias (__log1pl, log1pl)