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, see
  59     <http://www.gnu.org/licenses/>.  */
  60
  61 #include <math.h>
  62 #include <math_private.h>
  63
  64 /* log(1+x) = x - .5 x^2 + x^3 l(x)
  65    -.0078125 <= x <= +.0078125
  66    peak relative error 1.2e-37 */
  67 static const long double
  68 l3 =   3.333333333333333333333333333333336096926E-1L,
  69 l4 =  -2.499999999999999999999999999486853077002E-1L,
  70 l5 =   1.999999999999999999999999998515277861905E-1L,
  71 l6 =  -1.666666666666666666666798448356171665678E-1L,
  72 l7 =   1.428571428571428571428808945895490721564E-1L,
  73 l8 =  -1.249999999999999987884655626377588149000E-1L,
  74 l9 =   1.111111111111111093947834982832456459186E-1L,
  75 l10 = -1.000000000000532974938900317952530453248E-1L,
  76 l11 =  9.090909090915566247008015301349979892689E-2L,
  77 l12 = -8.333333211818065121250921925397567745734E-2L,
  78 l13 =  7.692307559897661630807048686258659316091E-2L,
  79 l14 = -7.144242754190814657241902218399056829264E-2L,
  80 l15 =  6.668057591071739754844678883223432347481E-2L;
  81
  82 /* Lookup table of ln(t) - (t-1)
  83     t = 0.5 + (k+26)/128)
  84     k = 0, ..., 91   */
  85 static const long double logtbl[92] = {
  86 -5.5345593589352099112142921677820359632418E-2L,
  87 -5.2108257402767124761784665198737642086148E-2L,
  88 -4.8991686870576856279407775480686721935120E-2L,
  89 -4.5993270766361228596215288742353061431071E-2L,
  90 -4.3110481649613269682442058976885699556950E-2L,
  91 -4.0340872319076331310838085093194799765520E-2L,
  92 -3.7682072451780927439219005993827431503510E-2L,
  93 -3.5131785416234343803903228503274262719586E-2L,
  94 -3.2687785249045246292687241862699949178831E-2L,
  95 -3.0347913785027239068190798397055267411813E-2L,
  96 -2.8110077931525797884641940838507561326298E-2L,
  97 -2.5972247078357715036426583294246819637618E-2L,
  98 -2.3932450635346084858612873953407168217307E-2L,
  99 -2.1988775689981395152022535153795155900240E-2L,
 100 -2.0139364778244501615441044267387667496733E-2L,
 101 -1.8382413762093794819267536615342902718324E-2L,
 102 -1.6716169807550022358923589720001638093023E-2L,
 103 -1.5138929457710992616226033183958974965355E-2L,
 104 -1.3649036795397472900424896523305726435029E-2L,
 105 -1.2244881690473465543308397998034325468152E-2L,
 106 -1.0924898127200937840689817557742469105693E-2L,
 107 -9.6875626072830301572839422532631079809328E-3L,
 108 -8.5313926245226231463436209313499745894157E-3L,
 109 -7.4549452072765973384933565912143044991706E-3L,
 110 -6.4568155251217050991200599386801665681310E-3L,
 111 -5.5356355563671005131126851708522185605193E-3L,
 112 -4.6900728132525199028885749289712348829878E-3L,
 113 -3.9188291218610470766469347968659624282519E-3L,
 114 -3.2206394539524058873423550293617843896540E-3L,
 115 -2.5942708080877805657374888909297113032132E-3L,
 116 -2.0385211375711716729239156839929281289086E-3L,
 117 -1.5522183228760777967376942769773768850872E-3L,
 118 -1.1342191863606077520036253234446621373191E-3L,
 119 -7.8340854719967065861624024730268350459991E-4L,
 120 -4.9869831458030115699628274852562992756174E-4L,
 121 -2.7902661731604211834685052867305795169688E-4L,
 122 -1.2335696813916860754951146082826952093496E-4L,
 123 -3.0677461025892873184042490943581654591817E-5L,
 124 #define ZERO logtbl[38]
 125  0.0000000000000000000000000000000000000000E0L,
 126 -3.0359557945051052537099938863236321874198E-5L,
 127 -1.2081346403474584914595395755316412213151E-4L,
 128 -2.7044071846562177120083903771008342059094E-4L,
 129 -4.7834133324631162897179240322783590830326E-4L,
 130 -7.4363569786340080624467487620270965403695E-4L,
 131 -1.0654639687057968333207323853366578860679E-3L,
 132 -1.4429854811877171341298062134712230604279E-3L,
 133 -1.8753781835651574193938679595797367137975E-3L,
 134 -2.3618380914922506054347222273705859653658E-3L,
 135 -2.9015787624124743013946600163375853631299E-3L,
 136 -3.4938307889254087318399313316921940859043E-3L,
 137 -4.1378413103128673800485306215154712148146E-3L,
 138 -4.8328735414488877044289435125365629849599E-3L,
 139 -5.5782063183564351739381962360253116934243E-3L,
 140 -6.3731336597098858051938306767880719015261E-3L,
 141 -7.2169643436165454612058905294782949315193E-3L,
 142 -8.1090214990427641365934846191367315083867E-3L,
 143 -9.0486422112807274112838713105168375482480E-3L,
 144 -1.0035177140880864314674126398350812606841E-2L,
 145 -1.1067990155502102718064936259435676477423E-2L,
 146 -1.2146457974158024928196575103115488672416E-2L,
 147 -1.3269969823361415906628825374158424754308E-2L,
 148 -1.4437927104692837124388550722759686270765E-2L,
 149 -1.5649743073340777659901053944852735064621E-2L,
 150 -1.6904842527181702880599758489058031645317E-2L,
 151 -1.8202661505988007336096407340750378994209E-2L,
 152 -1.9542647000370545390701192438691126552961E-2L,
 153 -2.0924256670080119637427928803038530924742E-2L,
 154 -2.2346958571309108496179613803760727786257E-2L,
 155 -2.3810230892650362330447187267648486279460E-2L,
 156 -2.5313561699385640380910474255652501521033E-2L,
 157 -2.6856448685790244233704909690165496625399E-2L,
 158 -2.8438398935154170008519274953860128449036E-2L,
 159 -3.0058928687233090922411781058956589863039E-2L,
 160 -3.1717563112854831855692484086486099896614E-2L,
 161 -3.3413836095418743219397234253475252001090E-2L,
 162 -3.5147290019036555862676702093393332533702E-2L,
 163 -3.6917475563073933027920505457688955423688E-2L,
 164 -3.8723951502862058660874073462456610731178E-2L,
 165 -4.0566284516358241168330505467000838017425E-2L,
 166 -4.2444048996543693813649967076598766917965E-2L,
 167 -4.4356826869355401653098777649745233339196E-2L,
 168 -4.6304207416957323121106944474331029996141E-2L,
 169 -4.8285787106164123613318093945035804818364E-2L,
 170 -5.0301169421838218987124461766244507342648E-2L,
 171 -5.2349964705088137924875459464622098310997E-2L,
 172 -5.4431789996103111613753440311680967840214E-2L,
 173 -5.6546268881465384189752786409400404404794E-2L,
 174 -5.8693031345788023909329239565012647817664E-2L,
 175 -6.0871713627532018185577188079210189048340E-2L,
 176 -6.3081958078862169742820420185833800925568E-2L,
 177 -6.5323413029406789694910800219643791556918E-2L,
 178 -6.7595732653791419081537811574227049288168E-2L
 179 };
 180
 181 /* ln(2) = ln2a + ln2b with extended precision. */
 182 static const long double
 183   ln2a = 6.93145751953125e-1L,
 184   ln2b = 1.4286068203094172321214581765680755001344E-6L;
 185
 186 static const long double
 187   ldbl_epsilon = 0x1p-106L;
 188
 189 long double
 190 __ieee754_logl(long double x)
 191 {
 192   long double z, y, w, t;
 193   unsigned int m;
 194   int k, e;
 195   double xhi;
 196   uint32_t hx, lx;
 197
 198   xhi = ldbl_high (x);
 199   EXTRACT_WORDS (hx, lx, xhi);
 200   m = hx;
 201
 202   /* Check for IEEE special cases.  */
 203   k = m & 0x7fffffff;
 204   /* log(0) = -infinity. */
 205   if ((k | lx) == 0)
 206     {
 207       return -0.5L / ZERO;
 208     }
 209   /* log ( x < 0 ) = NaN */
 210   if (m & 0x80000000)
 211     {
 212       return (x - x) / ZERO;
 213     }
 214   /* log (infinity or NaN) */
 215   if (k >= 0x7ff00000)
 216     {
 217       return x + x;
 218     }
 219
 220   /* On this interval the table is not used due to cancellation error.  */
 221   if ((x <= 1.0078125L) && (x >= 0.9921875L))
 222     {
 223       if (x == 1.0L)
 224         return 0.0L;
 225       z = x - 1.0L;
 226       k = 64;
 227       t = 1.0L;
 228       e = 0;
 229     }
 230   else
 231     {
 232       /* Extract exponent and reduce domain to 0.703125 <= u < 1.40625  */
 233       unsigned int w0;
 234       e = (int) (m >> 20) - (int) 0x3fe;
 235       if (e == -1022)
 236         {
 237           x *= 0x1p106L;
 238           xhi = ldbl_high (x);
 239           EXTRACT_WORDS (hx, lx, xhi);
 240           m = hx;
 241           e = (int) (m >> 20) - (int) 0x3fe - 106;
 242         }
 243       m &= 0xfffff;
 244       w0 = m | 0x3fe00000;
 245       m |= 0x100000;
 246       /* Find lookup table index k from high order bits of the significand. */
 247       if (m < 0x168000)
 248         {
 249           k = (m - 0xff000) >> 13;
 250           /* t is the argument 0.5 + (k+26)/128
 251              of the nearest item to u in the lookup table.  */
 252           INSERT_WORDS (xhi, 0x3ff00000 + (k << 13), 0);
 253           t = xhi;
 254           w0 += 0x100000;
 255           e -= 1;
 256           k += 64;
 257         }
 258       else
 259         {
 260           k = (m - 0xfe000) >> 14;
 261           INSERT_WORDS (xhi, 0x3fe00000 + (k << 14), 0);
 262           t = xhi;
 263         }
 264       x = __scalbnl (x, ((int) ((w0 - hx) * 2)) >> 21);
 265       /* log(u) = log( t u/t ) = log(t) + log(u/t)
 266          log(t) is tabulated in the lookup table.
 267          Express log(u/t) = log(1+z),  where z = u/t - 1 = (u-t)/t.
 268          cf. Cody & Waite. */
 269       z = (x - t) / t;
 270     }
 271   /* Series expansion of log(1+z).  */
 272   w = z * z;
 273   /* Avoid spurious underflows.  */
 274   if (__glibc_unlikely (fabsl (z) <= ldbl_epsilon))
 275     y = 0.0L;
 276   else
 277     {
 278       y = ((((((((((((l15 * z
 279                   + l14) * z
 280                  + l13) * z
 281                 + l12) * z
 282                + l11) * z
 283               + l10) * z
 284              + l9) * z
 285             + l8) * z
 286            + l7) * z
 287           + l6) * z
 288          + l5) * z
 289         + l4) * z
 290        + l3) * z * w;
 291       y -= 0.5 * w;
 292     }
 293   y += e * ln2b;  /* Base 2 exponent offset times ln(2).  */
 294   y += z;
 295   y += logtbl[k-26]; /* log(t) - (t-1) */
 296   y += (t - 1.0L);
 297   y += e * ln2a;
 298   return y;
 299 }
 300 strong_alias (__ieee754_logl, __logl_finite)