Back to index

glibc  2.9
e_log10l.c
Go to the documentation of this file.
00001 /*                                               log10l.c
00002  *
00003  *     Common logarithm, 128-bit long double precision
00004  *
00005  *
00006  *
00007  * SYNOPSIS:
00008  *
00009  * long double x, y, log10l();
00010  *
00011  * y = log10l( x );
00012  *
00013  *
00014  *
00015  * DESCRIPTION:
00016  *
00017  * Returns the base 10 logarithm of x.
00018  *
00019  * The argument is separated into its exponent and fractional
00020  * parts.  If the exponent is between -1 and +1, the logarithm
00021  * of the fraction is approximated by
00022  *
00023  *     log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
00024  *
00025  * Otherwise, setting  z = 2(x-1)/x+1),
00026  *
00027  *     log(x) = z + z^3 P(z)/Q(z).
00028  *
00029  *
00030  *
00031  * ACCURACY:
00032  *
00033  *                      Relative error:
00034  * arithmetic   domain     # trials      peak         rms
00035  *    IEEE      0.5, 2.0     30000      2.3e-34     4.9e-35
00036  *    IEEE     exp(+-10000)  30000      1.0e-34     4.1e-35
00037  *
00038  * In the tests over the interval exp(+-10000), the logarithms
00039  * of the random arguments were uniformly distributed over
00040  * [-10000, +10000].
00041  *
00042  */
00043 
00044 /*
00045    Cephes Math Library Release 2.2:  January, 1991
00046    Copyright 1984, 1991 by Stephen L. Moshier
00047    Adapted for glibc November, 2001
00048 
00049     This library is free software; you can redistribute it and/or
00050     modify it under the terms of the GNU Lesser General Public
00051     License as published by the Free Software Foundation; either
00052     version 2.1 of the License, or (at your option) any later version.
00053 
00054     This library is distributed in the hope that it will be useful,
00055     but WITHOUT ANY WARRANTY; without even the implied warranty of
00056     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00057     Lesser General Public License for more details.
00058 
00059     You should have received a copy of the GNU Lesser General Public
00060     License along with this library; if not, write to the Free Software
00061     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA 
00062 
00063  */
00064 
00065 #include "math.h"
00066 #include "math_private.h"
00067 
00068 /* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
00069  * 1/sqrt(2) <= x < sqrt(2)
00070  * Theoretical peak relative error = 5.3e-37,
00071  * relative peak error spread = 2.3e-14
00072  */
00073 static const long double P[13] =
00074 {
00075   1.313572404063446165910279910527789794488E4L,
00076   7.771154681358524243729929227226708890930E4L,
00077   2.014652742082537582487669938141683759923E5L,
00078   3.007007295140399532324943111654767187848E5L,
00079   2.854829159639697837788887080758954924001E5L,
00080   1.797628303815655343403735250238293741397E5L,
00081   7.594356839258970405033155585486712125861E4L,
00082   2.128857716871515081352991964243375186031E4L,
00083   3.824952356185897735160588078446136783779E3L,
00084   4.114517881637811823002128927449878962058E2L,
00085   2.321125933898420063925789532045674660756E1L,
00086   4.998469661968096229986658302195402690910E-1L,
00087   1.538612243596254322971797716843006400388E-6L
00088 };
00089 static const long double Q[12] =
00090 {
00091   3.940717212190338497730839731583397586124E4L,
00092   2.626900195321832660448791748036714883242E5L,
00093   7.777690340007566932935753241556479363645E5L,
00094   1.347518538384329112529391120390701166528E6L,
00095   1.514882452993549494932585972882995548426E6L,
00096   1.158019977462989115839826904108208787040E6L,
00097   6.132189329546557743179177159925690841200E5L,
00098   2.248234257620569139969141618556349415120E5L,
00099   5.605842085972455027590989944010492125825E4L,
00100   9.147150349299596453976674231612674085381E3L,
00101   9.104928120962988414618126155557301584078E2L,
00102   4.839208193348159620282142911143429644326E1L
00103 /* 1.000000000000000000000000000000000000000E0L, */
00104 };
00105 
00106 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
00107  * where z = 2(x-1)/(x+1)
00108  * 1/sqrt(2) <= x < sqrt(2)
00109  * Theoretical peak relative error = 1.1e-35,
00110  * relative peak error spread 1.1e-9
00111  */
00112 static const long double R[6] =
00113 {
00114   1.418134209872192732479751274970992665513E5L,
00115  -8.977257995689735303686582344659576526998E4L,
00116   2.048819892795278657810231591630928516206E4L,
00117  -2.024301798136027039250415126250455056397E3L,
00118   8.057002716646055371965756206836056074715E1L,
00119  -8.828896441624934385266096344596648080902E-1L
00120 };
00121 static const long double S[6] =
00122 {
00123   1.701761051846631278975701529965589676574E6L,
00124  -1.332535117259762928288745111081235577029E6L,
00125   4.001557694070773974936904547424676279307E5L,
00126  -5.748542087379434595104154610899551484314E4L,
00127   3.998526750980007367835804959888064681098E3L,
00128  -1.186359407982897997337150403816839480438E2L
00129 /* 1.000000000000000000000000000000000000000E0L, */
00130 };
00131 
00132 static const long double
00133 /* log10(2) */
00134 L102A = 0.3125L,
00135 L102B = -1.14700043360188047862611052755069732318101185E-2L,
00136 /* log10(e) */
00137 L10EA = 0.5L,
00138 L10EB = -6.570551809674817234887108108339491770560299E-2L,
00139 /* sqrt(2)/2 */
00140 SQRTH = 7.071067811865475244008443621048490392848359E-1L;
00141 
00142 
00143 
00144 /* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
00145 
00146 static long double
00147 neval (long double x, const long double *p, int n)
00148 {
00149   long double y;
00150 
00151   p += n;
00152   y = *p--;
00153   do
00154     {
00155       y = y * x + *p--;
00156     }
00157   while (--n > 0);
00158   return y;
00159 }
00160 
00161 
00162 /* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
00163 
00164 static long double
00165 deval (long double x, const long double *p, int n)
00166 {
00167   long double y;
00168 
00169   p += n;
00170   y = x + *p--;
00171   do
00172     {
00173       y = y * x + *p--;
00174     }
00175   while (--n > 0);
00176   return y;
00177 }
00178 
00179 
00180 
00181 long double
00182 __ieee754_log10l (x)
00183      long double x;
00184 {
00185   long double z;
00186   long double y;
00187   int e;
00188   int64_t hx, lx;
00189 
00190 /* Test for domain */
00191   GET_LDOUBLE_WORDS64 (hx, lx, x);
00192   if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
00193     return (-1.0L / (x - x));
00194   if (hx < 0)
00195     return (x - x) / (x - x);
00196   if (hx >= 0x7fff000000000000LL)
00197     return (x + x);
00198 
00199 /* separate mantissa from exponent */
00200 
00201 /* Note, frexp is used so that denormal numbers
00202  * will be handled properly.
00203  */
00204   x = __frexpl (x, &e);
00205 
00206 
00207 /* logarithm using log(x) = z + z**3 P(z)/Q(z),
00208  * where z = 2(x-1)/x+1)
00209  */
00210   if ((e > 2) || (e < -2))
00211     {
00212       if (x < SQRTH)
00213        {                    /* 2( 2x-1 )/( 2x+1 ) */
00214          e -= 1;
00215          z = x - 0.5L;
00216          y = 0.5L * z + 0.5L;
00217        }
00218       else
00219        {                    /*  2 (x-1)/(x+1)   */
00220          z = x - 0.5L;
00221          z -= 0.5L;
00222          y = 0.5L * x + 0.5L;
00223        }
00224       x = z / y;
00225       z = x * x;
00226       y = x * (z * neval (z, R, 5) / deval (z, S, 5));
00227       goto done;
00228     }
00229 
00230 
00231 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
00232 
00233   if (x < SQRTH)
00234     {
00235       e -= 1;
00236       x = 2.0 * x - 1.0L;   /*  2x - 1  */
00237     }
00238   else
00239     {
00240       x = x - 1.0L;
00241     }
00242   z = x * x;
00243   y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
00244   y = y - 0.5 * z;
00245 
00246 done:
00247 
00248   /* Multiply log of fraction by log10(e)
00249    * and base 2 exponent by log10(2).
00250    */
00251   z = y * L10EB;
00252   z += x * L10EB;
00253   z += e * L102B;
00254   z += y * L10EA;
00255   z += x * L10EA;
00256   z += e * L102A;
00257   return (z);
00258 }