Back to index

glibc  2.9
e_logl.c
Go to the documentation of this file.
00001 /*                                               logll.c
00002  *
00003  * Natural logarithm for 128-bit long double precision.
00004  *
00005  *
00006  *
00007  * SYNOPSIS:
00008  *
00009  * long double x, y, logl();
00010  *
00011  * y = logl( x );
00012  *
00013  *
00014  *
00015  * DESCRIPTION:
00016  *
00017  * Returns the base e (2.718...) logarithm of x.
00018  *
00019  * The argument is separated into its exponent and fractional
00020  * parts.  Use of a lookup table increases the speed of the routine.
00021  * The program uses logarithms tabulated at intervals of 1/128 to
00022  * cover the domain from approximately 0.7 to 1.4.
00023  *
00024  * On the interval [-1/128, +1/128] the logarithm of 1+x is approximated by
00025  *     log(1+x) = x - 0.5 x^2 + x^3 P(x) .
00026  *
00027  *
00028  *
00029  * ACCURACY:
00030  *
00031  *                      Relative error:
00032  * arithmetic   domain     # trials      peak         rms
00033  *    IEEE   0.875, 1.125   100000      1.2e-34    4.1e-35
00034  *    IEEE   0.125, 8       100000      1.2e-34    4.1e-35
00035  *
00036  *
00037  * WARNING:
00038  *
00039  * This program uses integer operations on bit fields of floating-point
00040  * numbers.  It does not work with data structures other than the
00041  * structure assumed.
00042  *
00043  */
00044 
00045 /* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov> 
00046 
00047     This library is free software; you can redistribute it and/or
00048     modify it under the terms of the GNU Lesser General Public
00049     License as published by the Free Software Foundation; either
00050     version 2.1 of the License, or (at your option) any later version.
00051 
00052     This library is distributed in the hope that it will be useful,
00053     but WITHOUT ANY WARRANTY; without even the implied warranty of
00054     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00055     Lesser General Public License for more details.
00056 
00057     You should have received a copy of the GNU Lesser General Public
00058     License along with this library; if not, write to the Free Software
00059     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
00060 
00061 #include "math_private.h"
00062 
00063 /* log(1+x) = x - .5 x^2 + x^3 l(x)
00064    -.0078125 <= x <= +.0078125
00065    peak relative error 1.2e-37 */
00066 static const long double
00067 l3 =   3.333333333333333333333333333333336096926E-1L,
00068 l4 =  -2.499999999999999999999999999486853077002E-1L,
00069 l5 =   1.999999999999999999999999998515277861905E-1L,
00070 l6 =  -1.666666666666666666666798448356171665678E-1L,
00071 l7 =   1.428571428571428571428808945895490721564E-1L,
00072 l8 =  -1.249999999999999987884655626377588149000E-1L,
00073 l9 =   1.111111111111111093947834982832456459186E-1L,
00074 l10 = -1.000000000000532974938900317952530453248E-1L,
00075 l11 =  9.090909090915566247008015301349979892689E-2L,
00076 l12 = -8.333333211818065121250921925397567745734E-2L,
00077 l13 =  7.692307559897661630807048686258659316091E-2L,
00078 l14 = -7.144242754190814657241902218399056829264E-2L,
00079 l15 =  6.668057591071739754844678883223432347481E-2L;
00080 
00081 /* Lookup table of ln(t) - (t-1)
00082     t = 0.5 + (k+26)/128)
00083     k = 0, ..., 91   */
00084 static const long double logtbl[92] = {
00085 -5.5345593589352099112142921677820359632418E-2L,
00086 -5.2108257402767124761784665198737642086148E-2L,
00087 -4.8991686870576856279407775480686721935120E-2L,
00088 -4.5993270766361228596215288742353061431071E-2L,
00089 -4.3110481649613269682442058976885699556950E-2L,
00090 -4.0340872319076331310838085093194799765520E-2L,
00091 -3.7682072451780927439219005993827431503510E-2L,
00092 -3.5131785416234343803903228503274262719586E-2L,
00093 -3.2687785249045246292687241862699949178831E-2L,
00094 -3.0347913785027239068190798397055267411813E-2L,
00095 -2.8110077931525797884641940838507561326298E-2L,
00096 -2.5972247078357715036426583294246819637618E-2L,
00097 -2.3932450635346084858612873953407168217307E-2L,
00098 -2.1988775689981395152022535153795155900240E-2L,
00099 -2.0139364778244501615441044267387667496733E-2L,
00100 -1.8382413762093794819267536615342902718324E-2L,
00101 -1.6716169807550022358923589720001638093023E-2L,
00102 -1.5138929457710992616226033183958974965355E-2L,
00103 -1.3649036795397472900424896523305726435029E-2L,
00104 -1.2244881690473465543308397998034325468152E-2L,
00105 -1.0924898127200937840689817557742469105693E-2L,
00106 -9.6875626072830301572839422532631079809328E-3L,
00107 -8.5313926245226231463436209313499745894157E-3L,
00108 -7.4549452072765973384933565912143044991706E-3L,
00109 -6.4568155251217050991200599386801665681310E-3L,
00110 -5.5356355563671005131126851708522185605193E-3L,
00111 -4.6900728132525199028885749289712348829878E-3L,
00112 -3.9188291218610470766469347968659624282519E-3L,
00113 -3.2206394539524058873423550293617843896540E-3L,
00114 -2.5942708080877805657374888909297113032132E-3L,
00115 -2.0385211375711716729239156839929281289086E-3L,
00116 -1.5522183228760777967376942769773768850872E-3L,
00117 -1.1342191863606077520036253234446621373191E-3L,
00118 -7.8340854719967065861624024730268350459991E-4L,
00119 -4.9869831458030115699628274852562992756174E-4L,
00120 -2.7902661731604211834685052867305795169688E-4L,
00121 -1.2335696813916860754951146082826952093496E-4L,
00122 -3.0677461025892873184042490943581654591817E-5L,
00123 #define ZERO logtbl[38]
00124  0.0000000000000000000000000000000000000000E0L,
00125 -3.0359557945051052537099938863236321874198E-5L,
00126 -1.2081346403474584914595395755316412213151E-4L,
00127 -2.7044071846562177120083903771008342059094E-4L,
00128 -4.7834133324631162897179240322783590830326E-4L,
00129 -7.4363569786340080624467487620270965403695E-4L,
00130 -1.0654639687057968333207323853366578860679E-3L,
00131 -1.4429854811877171341298062134712230604279E-3L,
00132 -1.8753781835651574193938679595797367137975E-3L,
00133 -2.3618380914922506054347222273705859653658E-3L,
00134 -2.9015787624124743013946600163375853631299E-3L,
00135 -3.4938307889254087318399313316921940859043E-3L,
00136 -4.1378413103128673800485306215154712148146E-3L,
00137 -4.8328735414488877044289435125365629849599E-3L,
00138 -5.5782063183564351739381962360253116934243E-3L,
00139 -6.3731336597098858051938306767880719015261E-3L,
00140 -7.2169643436165454612058905294782949315193E-3L,
00141 -8.1090214990427641365934846191367315083867E-3L,
00142 -9.0486422112807274112838713105168375482480E-3L,
00143 -1.0035177140880864314674126398350812606841E-2L,
00144 -1.1067990155502102718064936259435676477423E-2L,
00145 -1.2146457974158024928196575103115488672416E-2L,
00146 -1.3269969823361415906628825374158424754308E-2L,
00147 -1.4437927104692837124388550722759686270765E-2L,
00148 -1.5649743073340777659901053944852735064621E-2L,
00149 -1.6904842527181702880599758489058031645317E-2L,
00150 -1.8202661505988007336096407340750378994209E-2L,
00151 -1.9542647000370545390701192438691126552961E-2L,
00152 -2.0924256670080119637427928803038530924742E-2L,
00153 -2.2346958571309108496179613803760727786257E-2L,
00154 -2.3810230892650362330447187267648486279460E-2L,
00155 -2.5313561699385640380910474255652501521033E-2L,
00156 -2.6856448685790244233704909690165496625399E-2L,
00157 -2.8438398935154170008519274953860128449036E-2L,
00158 -3.0058928687233090922411781058956589863039E-2L,
00159 -3.1717563112854831855692484086486099896614E-2L,
00160 -3.3413836095418743219397234253475252001090E-2L,
00161 -3.5147290019036555862676702093393332533702E-2L,
00162 -3.6917475563073933027920505457688955423688E-2L,
00163 -3.8723951502862058660874073462456610731178E-2L,
00164 -4.0566284516358241168330505467000838017425E-2L,
00165 -4.2444048996543693813649967076598766917965E-2L,
00166 -4.4356826869355401653098777649745233339196E-2L,
00167 -4.6304207416957323121106944474331029996141E-2L,
00168 -4.8285787106164123613318093945035804818364E-2L,
00169 -5.0301169421838218987124461766244507342648E-2L,
00170 -5.2349964705088137924875459464622098310997E-2L,
00171 -5.4431789996103111613753440311680967840214E-2L,
00172 -5.6546268881465384189752786409400404404794E-2L,
00173 -5.8693031345788023909329239565012647817664E-2L,
00174 -6.0871713627532018185577188079210189048340E-2L,
00175 -6.3081958078862169742820420185833800925568E-2L,
00176 -6.5323413029406789694910800219643791556918E-2L,
00177 -6.7595732653791419081537811574227049288168E-2L
00178 };
00179 
00180 /* ln(2) = ln2a + ln2b with extended precision. */
00181 static const long double
00182   ln2a = 6.93145751953125e-1L,
00183   ln2b = 1.4286068203094172321214581765680755001344E-6L;
00184 
00185 long double
00186 __ieee754_logl(long double x)
00187 {
00188   long double z, y, w;
00189   ieee854_long_double_shape_type u, t;
00190   unsigned int m;
00191   int k, e;
00192 
00193   u.value = x;
00194   m = u.parts32.w0;
00195 
00196   /* Check for IEEE special cases.  */
00197   k = m & 0x7fffffff;
00198   /* log(0) = -infinity. */
00199   if ((k | u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
00200     {
00201       return -0.5L / ZERO;
00202     }
00203   /* log ( x < 0 ) = NaN */
00204   if (m & 0x80000000)
00205     {
00206       return (x - x) / ZERO;
00207     }
00208   /* log (infinity or NaN) */
00209   if (k >= 0x7fff0000)
00210     {
00211       return x + x;
00212     }
00213 
00214   /* Extract exponent and reduce domain to 0.703125 <= u < 1.40625  */
00215   e = (int) (m >> 16) - (int) 0x3ffe;
00216   m &= 0xffff;
00217   u.parts32.w0 = m | 0x3ffe0000;
00218   m |= 0x10000;
00219   /* Find lookup table index k from high order bits of the significand. */
00220   if (m < 0x16800)
00221     {
00222       k = (m - 0xff00) >> 9;
00223       /* t is the argument 0.5 + (k+26)/128
00224         of the nearest item to u in the lookup table.  */
00225       t.parts32.w0 = 0x3fff0000 + (k << 9);
00226       t.parts32.w1 = 0;
00227       t.parts32.w2 = 0;
00228       t.parts32.w3 = 0;
00229       u.parts32.w0 += 0x10000;
00230       e -= 1;
00231       k += 64;
00232     }
00233   else
00234     {
00235       k = (m - 0xfe00) >> 10;
00236       t.parts32.w0 = 0x3ffe0000 + (k << 10);
00237       t.parts32.w1 = 0;
00238       t.parts32.w2 = 0;
00239       t.parts32.w3 = 0;
00240     }
00241   /* On this interval the table is not used due to cancellation error.  */
00242   if ((x <= 1.0078125L) && (x >= 0.9921875L))
00243     {
00244       z = x - 1.0L;
00245       k = 64;
00246       t.value  = 1.0L;
00247       e = 0;
00248     }
00249   else
00250     {
00251       /* log(u) = log( t u/t ) = log(t) + log(u/t)
00252         log(t) is tabulated in the lookup table.
00253         Express log(u/t) = log(1+z),  where z = u/t - 1 = (u-t)/t.
00254          cf. Cody & Waite. */
00255       z = (u.value - t.value) / t.value;
00256     }
00257   /* Series expansion of log(1+z).  */
00258   w = z * z;
00259   y = ((((((((((((l15 * z
00260                 + l14) * z
00261                + l13) * z
00262               + l12) * z
00263               + l11) * z
00264              + l10) * z
00265             + l9) * z
00266            + l8) * z
00267           + l7) * z
00268          + l6) * z
00269         + l5) * z
00270        + l4) * z
00271        + l3) * z * w;
00272   y -= 0.5 * w;
00273   y += e * ln2b;  /* Base 2 exponent offset times ln(2).  */
00274   y += z;
00275   y += logtbl[k-26]; /* log(t) - (t-1) */
00276   y += (t.value - 1.0L);
00277   y += e * ln2a;
00278   return y;
00279 }