lib/logl.c

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