sysdeps/ieee754/ldbl-128ibm/e_logl.c

   1 /*                                                      logll.c
   2  *
   3  * Natural logarithm for 128-bit long double precision.
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * long double x, y, logl();
  10  *
  11  * y = logl( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns the base e (2.718...) logarithm of x.
  18  *
  19  * The argument is separated into its exponent and fractional
  20  * parts.  Use of a lookup table increases the speed of the routine.
  21  * The program uses logarithms tabulated at intervals of 1/128 to
  22  * cover the domain from approximately 0.7 to 1.4.
  23  *
  24  * On the interval [-1/128, +1/128] the logarithm of 1+x is approximated by
  25  *     log(1+x) = x - 0.5 x^2 + x^3 P(x) .
  26  *
  27  *
  28  *
  29  * ACCURACY:
  30  *
  31  *                      Relative error:
  32  * arithmetic   domain     # trials      peak         rms
  33  *    IEEE   0.875, 1.125   100000      1.2e-34    4.1e-35
  34  *    IEEE   0.125, 8       100000      1.2e-34    4.1e-35
  35  *
  36  *
  37  * WARNING:
  38  *
  39  * This program uses integer operations on bit fields of floating-point
  40  * numbers.  It does not work with data structures other than the
  41  * structure assumed.
  42  *
  43  */
  44
  45 /* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>
  46
  47     This library is free software; you can redistribute it and/or
  48     modify it under the terms of the GNU Lesser General Public
  49     License as published by the Free Software Foundation; either
  50     version 2.1 of the License, or (at your option) any later version.
  51
  52     This library is distributed in the hope that it will be useful,
  53     but WITHOUT ANY WARRANTY; without even the implied warranty of
  54     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  55     Lesser General Public License for more details.
  56
  57     You should have received a copy of the GNU Lesser General Public
  58     License along with this library; if not, write to the Free Software
  59     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
  60
  61 #include "math_private.h"
  62
  63 /* log(1+x) = x - .5 x^2 + x^3 l(x)
  64    -.0078125 <= x <= +.0078125
  65    peak relative error 1.2e-37 */
  66 static const long double
  67 l3 =   3.333333333333333333333333333333336096926E-1L,
  68 l4 =  -2.499999999999999999999999999486853077002E-1L,
  69 l5 =   1.999999999999999999999999998515277861905E-1L,
  70 l6 =  -1.666666666666666666666798448356171665678E-1L,
  71 l7 =   1.428571428571428571428808945895490721564E-1L,
  72 l8 =  -1.249999999999999987884655626377588149000E-1L,
  73 l9 =   1.111111111111111093947834982832456459186E-1L,
  74 l10 = -1.000000000000532974938900317952530453248E-1L,
  75 l11 =  9.090909090915566247008015301349979892689E-2L,
  76 l12 = -8.333333211818065121250921925397567745734E-2L,
  77 l13 =  7.692307559897661630807048686258659316091E-2L,
  78 l14 = -7.144242754190814657241902218399056829264E-2L,
  79 l15 =  6.668057591071739754844678883223432347481E-2L;
  80
  81 /* Lookup table of ln(t) - (t-1)
  82     t = 0.5 + (k+26)/128)
  83     k = 0, ..., 91   */
  84 static const long double logtbl[92] = {
  85 -5.5345593589352099112142921677820359632418E-2L,
  86 -5.2108257402767124761784665198737642086148E-2L,
  87 -4.8991686870576856279407775480686721935120E-2L,
  88 -4.5993270766361228596215288742353061431071E-2L,
  89 -4.3110481649613269682442058976885699556950E-2L,
  90 -4.0340872319076331310838085093194799765520E-2L,
  91 -3.7682072451780927439219005993827431503510E-2L,
  92 -3.5131785416234343803903228503274262719586E-2L,
  93 -3.2687785249045246292687241862699949178831E-2L,
  94 -3.0347913785027239068190798397055267411813E-2L,
  95 -2.8110077931525797884641940838507561326298E-2L,
  96 -2.5972247078357715036426583294246819637618E-2L,
  97 -2.3932450635346084858612873953407168217307E-2L,
  98 -2.1988775689981395152022535153795155900240E-2L,
  99 -2.0139364778244501615441044267387667496733E-2L,
 100 -1.8382413762093794819267536615342902718324E-2L,
 101 -1.6716169807550022358923589720001638093023E-2L,
 102 -1.5138929457710992616226033183958974965355E-2L,
 103 -1.3649036795397472900424896523305726435029E-2L,
 104 -1.2244881690473465543308397998034325468152E-2L,
 105 -1.0924898127200937840689817557742469105693E-2L,
 106 -9.6875626072830301572839422532631079809328E-3L,
 107 -8.5313926245226231463436209313499745894157E-3L,
 108 -7.4549452072765973384933565912143044991706E-3L,
 109 -6.4568155251217050991200599386801665681310E-3L,
 110 -5.5356355563671005131126851708522185605193E-3L,
 111 -4.6900728132525199028885749289712348829878E-3L,
 112 -3.9188291218610470766469347968659624282519E-3L,
 113 -3.2206394539524058873423550293617843896540E-3L,
 114 -2.5942708080877805657374888909297113032132E-3L,
 115 -2.0385211375711716729239156839929281289086E-3L,
 116 -1.5522183228760777967376942769773768850872E-3L,
 117 -1.1342191863606077520036253234446621373191E-3L,
 118 -7.8340854719967065861624024730268350459991E-4L,
 119 -4.9869831458030115699628274852562992756174E-4L,
 120 -2.7902661731604211834685052867305795169688E-4L,
 121 -1.2335696813916860754951146082826952093496E-4L,
 122 -3.0677461025892873184042490943581654591817E-5L,
 123 #define ZERO logtbl[38]
 124  0.0000000000000000000000000000000000000000E0L,
 125 -3.0359557945051052537099938863236321874198E-5L,
 126 -1.2081346403474584914595395755316412213151E-4L,
 127 -2.7044071846562177120083903771008342059094E-4L,
 128 -4.7834133324631162897179240322783590830326E-4L,
 129 -7.4363569786340080624467487620270965403695E-4L,
 130 -1.0654639687057968333207323853366578860679E-3L,
 131 -1.4429854811877171341298062134712230604279E-3L,
 132 -1.8753781835651574193938679595797367137975E-3L,
 133 -2.3618380914922506054347222273705859653658E-3L,
 134 -2.9015787624124743013946600163375853631299E-3L,
 135 -3.4938307889254087318399313316921940859043E-3L,
 136 -4.1378413103128673800485306215154712148146E-3L,
 137 -4.8328735414488877044289435125365629849599E-3L,
 138 -5.5782063183564351739381962360253116934243E-3L,
 139 -6.3731336597098858051938306767880719015261E-3L,
 140 -7.2169643436165454612058905294782949315193E-3L,
 141 -8.1090214990427641365934846191367315083867E-3L,
 142 -9.0486422112807274112838713105168375482480E-3L,
 143 -1.0035177140880864314674126398350812606841E-2L,
 144 -1.1067990155502102718064936259435676477423E-2L,
 145 -1.2146457974158024928196575103115488672416E-2L,
 146 -1.3269969823361415906628825374158424754308E-2L,
 147 -1.4437927104692837124388550722759686270765E-2L,
 148 -1.5649743073340777659901053944852735064621E-2L,
 149 -1.6904842527181702880599758489058031645317E-2L,
 150 -1.8202661505988007336096407340750378994209E-2L,
 151 -1.9542647000370545390701192438691126552961E-2L,
 152 -2.0924256670080119637427928803038530924742E-2L,
 153 -2.2346958571309108496179613803760727786257E-2L,
 154 -2.3810230892650362330447187267648486279460E-2L,
 155 -2.5313561699385640380910474255652501521033E-2L,
 156 -2.6856448685790244233704909690165496625399E-2L,
 157 -2.8438398935154170008519274953860128449036E-2L,
 158 -3.0058928687233090922411781058956589863039E-2L,
 159 -3.1717563112854831855692484086486099896614E-2L,
 160 -3.3413836095418743219397234253475252001090E-2L,
 161 -3.5147290019036555862676702093393332533702E-2L,
 162 -3.6917475563073933027920505457688955423688E-2L,
 163 -3.8723951502862058660874073462456610731178E-2L,
 164 -4.0566284516358241168330505467000838017425E-2L,
 165 -4.2444048996543693813649967076598766917965E-2L,
 166 -4.4356826869355401653098777649745233339196E-2L,
 167 -4.6304207416957323121106944474331029996141E-2L,
 168 -4.8285787106164123613318093945035804818364E-2L,
 169 -5.0301169421838218987124461766244507342648E-2L,
 170 -5.2349964705088137924875459464622098310997E-2L,
 171 -5.4431789996103111613753440311680967840214E-2L,
 172 -5.6546268881465384189752786409400404404794E-2L,
 173 -5.8693031345788023909329239565012647817664E-2L,
 174 -6.0871713627532018185577188079210189048340E-2L,
 175 -6.3081958078862169742820420185833800925568E-2L,
 176 -6.5323413029406789694910800219643791556918E-2L,
 177 -6.7595732653791419081537811574227049288168E-2L
 178 };
 179
 180 /* ln(2) = ln2a + ln2b with extended precision. */
 181 static const long double
 182   ln2a = 6.93145751953125e-1L,
 183   ln2b = 1.4286068203094172321214581765680755001344E-6L;
 184
 185 long double
 186 __ieee754_logl(long double x)
 187 {
 188   long double z, y, w;
 189   ieee854_long_double_shape_type u, t;
 190   unsigned int m;
 191   int k, e;
 192
 193   u.value = x;
 194   m = u.parts32.w0;
 195
 196   /* Check for IEEE special cases.  */
 197   k = m & 0x7fffffff;
 198   /* log(0) = -infinity. */
 199   if ((k | u.parts32.w1 | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3) == 0)
 200     {
 201       return -0.5L / ZERO;
 202     }
 203   /* log ( x < 0 ) = NaN */
 204   if (m & 0x80000000)
 205     {
 206       return (x - x) / ZERO;
 207     }
 208   /* log (infinity or NaN) */
 209   if (k >= 0x7ff00000)
 210     {
 211       return x + x;
 212     }
 213
 214   /* On this interval the table is not used due to cancellation error.  */
 215   if ((x <= 1.0078125L) && (x >= 0.9921875L))
 216     {
 217       z = x - 1.0L;
 218       k = 64;
 219       t.value  = 1.0L;
 220       e = 0;
 221     }
 222   else
 223     {
 224       /* Extract exponent and reduce domain to 0.703125 <= u < 1.40625  */
 225       unsigned int w0;
 226       e = (int) (m >> 20) - (int) 0x3fe;
 227       m &= 0xfffff;
 228       w0 = m | 0x3fe00000;
 229       m |= 0x100000;
 230       /* Find lookup table index k from high order bits of the significand. */
 231       if (m < 0x168000)
 232         {
 233           k = (m - 0xff000) >> 13;
 234           /* t is the argument 0.5 + (k+26)/128
 235              of the nearest item to u in the lookup table.  */
 236           t.parts32.w0 = 0x3ff00000 + (k << 13);
 237           t.parts32.w1 = 0;
 238           t.parts32.w2 = 0;
 239           t.parts32.w3 = 0;
 240           w0 += 0x100000;
 241           e -= 1;
 242           k += 64;
 243         }
 244       else
 245         {
 246           k = (m - 0xfe000) >> 14;
 247           t.parts32.w0 = 0x3fe00000 + (k << 14);
 248           t.parts32.w1 = 0;
 249           t.parts32.w2 = 0;
 250           t.parts32.w3 = 0;
 251         }
 252       u.value = __scalbnl (u.value, ((int) ((w0 - u.parts32.w0) * 2)) >> 21);
 253       /* log(u) = log( t u/t ) = log(t) + log(u/t)
 254          log(t) is tabulated in the lookup table.
 255          Express log(u/t) = log(1+z),  where z = u/t - 1 = (u-t)/t.
 256          cf. Cody & Waite. */
 257       z = (u.value - t.value) / t.value;
 258     }
 259   /* Series expansion of log(1+z).  */
 260   w = z * z;
 261   y = ((((((((((((l15 * z
 262                   + l14) * z
 263                  + l13) * z
 264                 + l12) * z
 265                + l11) * z
 266               + l10) * z
 267              + l9) * z
 268             + l8) * z
 269            + l7) * z
 270           + l6) * z
 271          + l5) * z
 272         + l4) * z
 273        + l3) * z * w;
 274   y -= 0.5 * w;
 275   y += e * ln2b;  /* Base 2 exponent offset times ln(2).  */
 276   y += z;
 277   y += logtbl[k-26]; /* log(t) - (t-1) */
 278   y += (t.value - 1.0L);
 279   y += e * ln2a;
 280   return y;
 281 }