libquadmath/math/logq.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, logq();
  10  *
  11  * y = logq( 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 "quadmath-imp.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 __float128
  67 l3 =   3.333333333333333333333333333333336096926E-1Q,
  68 l4 =  -2.499999999999999999999999999486853077002E-1Q,
  69 l5 =   1.999999999999999999999999998515277861905E-1Q,
  70 l6 =  -1.666666666666666666666798448356171665678E-1Q,
  71 l7 =   1.428571428571428571428808945895490721564E-1Q,
  72 l8 =  -1.249999999999999987884655626377588149000E-1Q,
  73 l9 =   1.111111111111111093947834982832456459186E-1Q,
  74 l10 = -1.000000000000532974938900317952530453248E-1Q,
  75 l11 =  9.090909090915566247008015301349979892689E-2Q,
  76 l12 = -8.333333211818065121250921925397567745734E-2Q,
  77 l13 =  7.692307559897661630807048686258659316091E-2Q,
  78 l14 = -7.144242754190814657241902218399056829264E-2Q,
  79 l15 =  6.668057591071739754844678883223432347481E-2Q;
  80
  81 /* Lookup table of ln(t) - (t-1)
  82     t = 0.5 + (k+26)/128)
  83     k = 0, ..., 91   */
  84 static const __float128 logtbl[92] = {
  85 -5.5345593589352099112142921677820359632418E-2Q,
  86 -5.2108257402767124761784665198737642086148E-2Q,
  87 -4.8991686870576856279407775480686721935120E-2Q,
  88 -4.5993270766361228596215288742353061431071E-2Q,
  89 -4.3110481649613269682442058976885699556950E-2Q,
  90 -4.0340872319076331310838085093194799765520E-2Q,
  91 -3.7682072451780927439219005993827431503510E-2Q,
  92 -3.5131785416234343803903228503274262719586E-2Q,
  93 -3.2687785249045246292687241862699949178831E-2Q,
  94 -3.0347913785027239068190798397055267411813E-2Q,
  95 -2.8110077931525797884641940838507561326298E-2Q,
  96 -2.5972247078357715036426583294246819637618E-2Q,
  97 -2.3932450635346084858612873953407168217307E-2Q,
  98 -2.1988775689981395152022535153795155900240E-2Q,
  99 -2.0139364778244501615441044267387667496733E-2Q,
 100 -1.8382413762093794819267536615342902718324E-2Q,
 101 -1.6716169807550022358923589720001638093023E-2Q,
 102 -1.5138929457710992616226033183958974965355E-2Q,
 103 -1.3649036795397472900424896523305726435029E-2Q,
 104 -1.2244881690473465543308397998034325468152E-2Q,
 105 -1.0924898127200937840689817557742469105693E-2Q,
 106 -9.6875626072830301572839422532631079809328E-3Q,
 107 -8.5313926245226231463436209313499745894157E-3Q,
 108 -7.4549452072765973384933565912143044991706E-3Q,
 109 -6.4568155251217050991200599386801665681310E-3Q,
 110 -5.5356355563671005131126851708522185605193E-3Q,
 111 -4.6900728132525199028885749289712348829878E-3Q,
 112 -3.9188291218610470766469347968659624282519E-3Q,
 113 -3.2206394539524058873423550293617843896540E-3Q,
 114 -2.5942708080877805657374888909297113032132E-3Q,
 115 -2.0385211375711716729239156839929281289086E-3Q,
 116 -1.5522183228760777967376942769773768850872E-3Q,
 117 -1.1342191863606077520036253234446621373191E-3Q,
 118 -7.8340854719967065861624024730268350459991E-4Q,
 119 -4.9869831458030115699628274852562992756174E-4Q,
 120 -2.7902661731604211834685052867305795169688E-4Q,
 121 -1.2335696813916860754951146082826952093496E-4Q,
 122 -3.0677461025892873184042490943581654591817E-5Q,
 123 #define ZERO logtbl[38]
 124  0.0000000000000000000000000000000000000000E0Q,
 125 -3.0359557945051052537099938863236321874198E-5Q,
 126 -1.2081346403474584914595395755316412213151E-4Q,
 127 -2.7044071846562177120083903771008342059094E-4Q,
 128 -4.7834133324631162897179240322783590830326E-4Q,
 129 -7.4363569786340080624467487620270965403695E-4Q,
 130 -1.0654639687057968333207323853366578860679E-3Q,
 131 -1.4429854811877171341298062134712230604279E-3Q,
 132 -1.8753781835651574193938679595797367137975E-3Q,
 133 -2.3618380914922506054347222273705859653658E-3Q,
 134 -2.9015787624124743013946600163375853631299E-3Q,
 135 -3.4938307889254087318399313316921940859043E-3Q,
 136 -4.1378413103128673800485306215154712148146E-3Q,
 137 -4.8328735414488877044289435125365629849599E-3Q,
 138 -5.5782063183564351739381962360253116934243E-3Q,
 139 -6.3731336597098858051938306767880719015261E-3Q,
 140 -7.2169643436165454612058905294782949315193E-3Q,
 141 -8.1090214990427641365934846191367315083867E-3Q,
 142 -9.0486422112807274112838713105168375482480E-3Q,
 143 -1.0035177140880864314674126398350812606841E-2Q,
 144 -1.1067990155502102718064936259435676477423E-2Q,
 145 -1.2146457974158024928196575103115488672416E-2Q,
 146 -1.3269969823361415906628825374158424754308E-2Q,
 147 -1.4437927104692837124388550722759686270765E-2Q,
 148 -1.5649743073340777659901053944852735064621E-2Q,
 149 -1.6904842527181702880599758489058031645317E-2Q,
 150 -1.8202661505988007336096407340750378994209E-2Q,
 151 -1.9542647000370545390701192438691126552961E-2Q,
 152 -2.0924256670080119637427928803038530924742E-2Q,
 153 -2.2346958571309108496179613803760727786257E-2Q,
 154 -2.3810230892650362330447187267648486279460E-2Q,
 155 -2.5313561699385640380910474255652501521033E-2Q,
 156 -2.6856448685790244233704909690165496625399E-2Q,
 157 -2.8438398935154170008519274953860128449036E-2Q,
 158 -3.0058928687233090922411781058956589863039E-2Q,
 159 -3.1717563112854831855692484086486099896614E-2Q,
 160 -3.3413836095418743219397234253475252001090E-2Q,
 161 -3.5147290019036555862676702093393332533702E-2Q,
 162 -3.6917475563073933027920505457688955423688E-2Q,
 163 -3.8723951502862058660874073462456610731178E-2Q,
 164 -4.0566284516358241168330505467000838017425E-2Q,
 165 -4.2444048996543693813649967076598766917965E-2Q,
 166 -4.4356826869355401653098777649745233339196E-2Q,
 167 -4.6304207416957323121106944474331029996141E-2Q,
 168 -4.8285787106164123613318093945035804818364E-2Q,
 169 -5.0301169421838218987124461766244507342648E-2Q,
 170 -5.2349964705088137924875459464622098310997E-2Q,
 171 -5.4431789996103111613753440311680967840214E-2Q,
 172 -5.6546268881465384189752786409400404404794E-2Q,
 173 -5.8693031345788023909329239565012647817664E-2Q,
 174 -6.0871713627532018185577188079210189048340E-2Q,
 175 -6.3081958078862169742820420185833800925568E-2Q,
 176 -6.5323413029406789694910800219643791556918E-2Q,
 177 -6.7595732653791419081537811574227049288168E-2Q
 178 };
 179
 180 /* ln(2) = ln2a + ln2b with extended precision. */
 181 static const __float128
 182   ln2a = 6.93145751953125e-1Q,
 183   ln2b = 1.4286068203094172321214581765680755001344E-6Q;
 184
 185 __float128
 186 logq(__float128 x)
 187 {
 188   __float128 z, y, w;
 189   ieee854_float128 u, t;
 190   unsigned int m;
 191   int k, e;
 192
 193   u.value = x;
 194   m = u.words32.w0;
 195
 196   /* Check for IEEE special cases.  */
 197   k = m & 0x7fffffff;
 198   /* log(0) = -infinity. */
 199   if ((k | u.words32.w1 | u.words32.w2 | u.words32.w3) == 0)
 200     {
 201       return -0.5Q / 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 >= 0x7fff0000)
 210     {
 211       return x + x;
 212     }
 213
 214   /* Extract exponent and reduce domain to 0.703125 <= u < 1.40625  */
 215   u.value = frexpq (x, &e);
 216   m = u.words32.w0 & 0xffff;
 217   m |= 0x10000;
 218   /* Find lookup table index k from high order bits of the significand. */
 219   if (m < 0x16800)
 220     {
 221       k = (m - 0xff00) >> 9;
 222       /* t is the argument 0.5 + (k+26)/128
 223          of the nearest item to u in the lookup table.  */
 224       t.words32.w0 = 0x3fff0000 + (k << 9);
 225       t.words32.w1 = 0;
 226       t.words32.w2 = 0;
 227       t.words32.w3 = 0;
 228       u.words32.w0 += 0x10000;
 229       e -= 1;
 230       k += 64;
 231     }
 232   else
 233     {
 234       k = (m - 0xfe00) >> 10;
 235       t.words32.w0 = 0x3ffe0000 + (k << 10);
 236       t.words32.w1 = 0;
 237       t.words32.w2 = 0;
 238       t.words32.w3 = 0;
 239     }
 240   /* On this interval the table is not used due to cancellation error.  */
 241   if ((x <= 1.0078125Q) && (x >= 0.9921875Q))
 242     {
 243       if (x == 1)
 244         return 0;
 245       z = x - 1;
 246       k = 64;
 247       t.value  = 1;
 248       e = 0;
 249     }
 250   else
 251     {
 252       /* log(u) = log( t u/t ) = log(t) + log(u/t)
 253          log(t) is tabulated in the lookup table.
 254          Express log(u/t) = log(1+z),  where z = u/t - 1 = (u-t)/t.
 255          cf. Cody & Waite. */
 256       z = (u.value - t.value) / t.value;
 257     }
 258   /* Series expansion of log(1+z).  */
 259   w = z * z;
 260   y = ((((((((((((l15 * z
 261                   + l14) * z
 262                  + l13) * z
 263                 + l12) * z
 264                + l11) * z
 265               + l10) * z
 266              + l9) * z
 267             + l8) * z
 268            + l7) * z
 269           + l6) * z
 270          + l5) * z
 271         + l4) * z
 272        + l3) * z * w;
 273   y -= 0.5 * w;
 274   y += e * ln2b;  /* Base 2 exponent offset times ln(2).  */
 275   y += z;
 276   y += logtbl[k-26]; /* log(t) - (t-1) */
 277   y += (t.value - 1);
 278   y += e * ln2a;
 279   return y;
 280 }