Back to index

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