1 /* $OpenBSD: e_logl.c,v 1.1 2011/07/06 00:02:42 martynas Exp $ */ 2 3 /* 4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> 5 * 6 * Permission to use, copy, modify, and distribute this software for any 7 * purpose with or without fee is hereby granted, provided that the above 8 * copyright notice and this permission notice appear in all copies. 9 * 10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES 11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF 12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR 13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN 15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 17 */ 18 19 /* logl.c 20 * 21 * Natural logarithm for 128-bit long double precision. 22 * 23 * 24 * 25 * SYNOPSIS: 26 * 27 * long double x, y, logl(); 28 * 29 * y = logl( x ); 30 * 31 * 32 * 33 * DESCRIPTION: 34 * 35 * Returns the base e (2.718...) logarithm of x. 36 * 37 * The argument is separated into its exponent and fractional 38 * parts. Use of a lookup table increases the speed of the routine. 39 * The program uses logarithms tabulated at intervals of 1/128 to 40 * cover the domain from approximately 0.7 to 1.4. 41 * 42 * On the interval [-1/128, +1/128] the logarithm of 1+x is approximated by 43 * log(1+x) = x - 0.5 x^2 + x^3 P(x) . 44 * 45 * 46 * 47 * ACCURACY: 48 * 49 * Relative error: 50 * arithmetic domain # trials peak rms 51 * IEEE 0.875, 1.125 100000 1.2e-34 4.1e-35 52 * IEEE 0.125, 8 100000 1.2e-34 4.1e-35 53 * 54 * 55 * WARNING: 56 * 57 * This program uses integer operations on bit fields of floating-point 58 * numbers. It does not work with data structures other than the 59 * structure assumed. 60 * 61 */ 62 63 #include <openlibm.h> 64 65 #include "math_private.h" 66 67 /* log(1+x) = x - .5 x^2 + x^3 l(x) 68 -.0078125 <= x <= +.0078125 69 peak relative error 1.2e-37 */ 70 static const long double 71 l3 = 3.333333333333333333333333333333336096926E-1L, 72 l4 = -2.499999999999999999999999999486853077002E-1L, 73 l5 = 1.999999999999999999999999998515277861905E-1L, 74 l6 = -1.666666666666666666666798448356171665678E-1L, 75 l7 = 1.428571428571428571428808945895490721564E-1L, 76 l8 = -1.249999999999999987884655626377588149000E-1L, 77 l9 = 1.111111111111111093947834982832456459186E-1L, 78 l10 = -1.000000000000532974938900317952530453248E-1L, 79 l11 = 9.090909090915566247008015301349979892689E-2L, 80 l12 = -8.333333211818065121250921925397567745734E-2L, 81 l13 = 7.692307559897661630807048686258659316091E-2L, 82 l14 = -7.144242754190814657241902218399056829264E-2L, 83 l15 = 6.668057591071739754844678883223432347481E-2L; 84 85 /* Lookup table of ln(t) - (t-1) 86 t = 0.5 + (k+26)/128) 87 k = 0, ..., 91 */ 88 static const long double logtbl[92] = { 89 -5.5345593589352099112142921677820359632418E-2L, 90 -5.2108257402767124761784665198737642086148E-2L, 91 -4.8991686870576856279407775480686721935120E-2L, 92 -4.5993270766361228596215288742353061431071E-2L, 93 -4.3110481649613269682442058976885699556950E-2L, 94 -4.0340872319076331310838085093194799765520E-2L, 95 -3.7682072451780927439219005993827431503510E-2L, 96 -3.5131785416234343803903228503274262719586E-2L, 97 -3.2687785249045246292687241862699949178831E-2L, 98 -3.0347913785027239068190798397055267411813E-2L, 99 -2.8110077931525797884641940838507561326298E-2L, 100 -2.5972247078357715036426583294246819637618E-2L, 101 -2.3932450635346084858612873953407168217307E-2L, 102 -2.1988775689981395152022535153795155900240E-2L, 103 -2.0139364778244501615441044267387667496733E-2L, 104 -1.8382413762093794819267536615342902718324E-2L, 105 -1.6716169807550022358923589720001638093023E-2L, 106 -1.5138929457710992616226033183958974965355E-2L, 107 -1.3649036795397472900424896523305726435029E-2L, 108 -1.2244881690473465543308397998034325468152E-2L, 109 -1.0924898127200937840689817557742469105693E-2L, 110 -9.6875626072830301572839422532631079809328E-3L, 111 -8.5313926245226231463436209313499745894157E-3L, 112 -7.4549452072765973384933565912143044991706E-3L, 113 -6.4568155251217050991200599386801665681310E-3L, 114 -5.5356355563671005131126851708522185605193E-3L, 115 -4.6900728132525199028885749289712348829878E-3L, 116 -3.9188291218610470766469347968659624282519E-3L, 117 -3.2206394539524058873423550293617843896540E-3L, 118 -2.5942708080877805657374888909297113032132E-3L, 119 -2.0385211375711716729239156839929281289086E-3L, 120 -1.5522183228760777967376942769773768850872E-3L, 121 -1.1342191863606077520036253234446621373191E-3L, 122 -7.8340854719967065861624024730268350459991E-4L, 123 -4.9869831458030115699628274852562992756174E-4L, 124 -2.7902661731604211834685052867305795169688E-4L, 125 -1.2335696813916860754951146082826952093496E-4L, 126 -3.0677461025892873184042490943581654591817E-5L, 127 #define ZERO logtbl[38] 128 0.0000000000000000000000000000000000000000E0L, 129 -3.0359557945051052537099938863236321874198E-5L, 130 -1.2081346403474584914595395755316412213151E-4L, 131 -2.7044071846562177120083903771008342059094E-4L, 132 -4.7834133324631162897179240322783590830326E-4L, 133 -7.4363569786340080624467487620270965403695E-4L, 134 -1.0654639687057968333207323853366578860679E-3L, 135 -1.4429854811877171341298062134712230604279E-3L, 136 -1.8753781835651574193938679595797367137975E-3L, 137 -2.3618380914922506054347222273705859653658E-3L, 138 -2.9015787624124743013946600163375853631299E-3L, 139 -3.4938307889254087318399313316921940859043E-3L, 140 -4.1378413103128673800485306215154712148146E-3L, 141 -4.8328735414488877044289435125365629849599E-3L, 142 -5.5782063183564351739381962360253116934243E-3L, 143 -6.3731336597098858051938306767880719015261E-3L, 144 -7.2169643436165454612058905294782949315193E-3L, 145 -8.1090214990427641365934846191367315083867E-3L, 146 -9.0486422112807274112838713105168375482480E-3L, 147 -1.0035177140880864314674126398350812606841E-2L, 148 -1.1067990155502102718064936259435676477423E-2L, 149 -1.2146457974158024928196575103115488672416E-2L, 150 -1.3269969823361415906628825374158424754308E-2L, 151 -1.4437927104692837124388550722759686270765E-2L, 152 -1.5649743073340777659901053944852735064621E-2L, 153 -1.6904842527181702880599758489058031645317E-2L, 154 -1.8202661505988007336096407340750378994209E-2L, 155 -1.9542647000370545390701192438691126552961E-2L, 156 -2.0924256670080119637427928803038530924742E-2L, 157 -2.2346958571309108496179613803760727786257E-2L, 158 -2.3810230892650362330447187267648486279460E-2L, 159 -2.5313561699385640380910474255652501521033E-2L, 160 -2.6856448685790244233704909690165496625399E-2L, 161 -2.8438398935154170008519274953860128449036E-2L, 162 -3.0058928687233090922411781058956589863039E-2L, 163 -3.1717563112854831855692484086486099896614E-2L, 164 -3.3413836095418743219397234253475252001090E-2L, 165 -3.5147290019036555862676702093393332533702E-2L, 166 -3.6917475563073933027920505457688955423688E-2L, 167 -3.8723951502862058660874073462456610731178E-2L, 168 -4.0566284516358241168330505467000838017425E-2L, 169 -4.2444048996543693813649967076598766917965E-2L, 170 -4.4356826869355401653098777649745233339196E-2L, 171 -4.6304207416957323121106944474331029996141E-2L, 172 -4.8285787106164123613318093945035804818364E-2L, 173 -5.0301169421838218987124461766244507342648E-2L, 174 -5.2349964705088137924875459464622098310997E-2L, 175 -5.4431789996103111613753440311680967840214E-2L, 176 -5.6546268881465384189752786409400404404794E-2L, 177 -5.8693031345788023909329239565012647817664E-2L, 178 -6.0871713627532018185577188079210189048340E-2L, 179 -6.3081958078862169742820420185833800925568E-2L, 180 -6.5323413029406789694910800219643791556918E-2L, 181 -6.7595732653791419081537811574227049288168E-2L 182 }; 183 184 /* ln(2) = ln2a + ln2b with extended precision. */ 185 static const long double 186 ln2a = 6.93145751953125e-1L, 187 ln2b = 1.4286068203094172321214581765680755001344E-6L; 188 189 long double 190 logl(long double x) 191 { 192 long double z, y, w; 193 ieee_quad_shape_type u, t; 194 unsigned int m; 195 int k, e; 196 197 u.value = x; 198 m = u.parts32.mswhi; 199 200 /* Check for IEEE special cases. */ 201 k = m & 0x7fffffff; 202 /* log(0) = -infinity. */ 203 if ((k | u.parts32.mswlo | u.parts32.lswhi | u.parts32.lswlo) == 0) 204 { 205 return -0.5L / ZERO; 206 } 207 /* log ( x < 0 ) = NaN */ 208 if (m & 0x80000000) 209 { 210 return (x - x) / ZERO; 211 } 212 /* log (infinity or NaN) */ 213 if (k >= 0x7fff0000) 214 { 215 return x + x; 216 } 217 218 /* Extract exponent and reduce domain to 0.703125 <= u < 1.40625 */ 219 e = (int) (m >> 16) - (int) 0x3ffe; 220 m &= 0xffff; 221 u.parts32.mswhi = m | 0x3ffe0000; 222 m |= 0x10000; 223 /* Find lookup table index k from high order bits of the significand. */ 224 if (m < 0x16800) 225 { 226 k = (m - 0xff00) >> 9; 227 /* t is the argument 0.5 + (k+26)/128 228 of the nearest item to u in the lookup table. */ 229 t.parts32.mswhi = 0x3fff0000 + (k << 9); 230 t.parts32.mswlo = 0; 231 t.parts32.lswhi = 0; 232 t.parts32.lswlo = 0; 233 u.parts32.mswhi += 0x10000; 234 e -= 1; 235 k += 64; 236 } 237 else 238 { 239 k = (m - 0xfe00) >> 10; 240 t.parts32.mswhi = 0x3ffe0000 + (k << 10); 241 t.parts32.mswlo = 0; 242 t.parts32.lswhi = 0; 243 t.parts32.lswlo = 0; 244 } 245 /* On this interval the table is not used due to cancellation error. */ 246 if ((x <= 1.0078125L) && (x >= 0.9921875L)) 247 { 248 z = x - 1.0L; 249 k = 64; 250 t.value = 1.0L; 251 e = 0; 252 } 253 else 254 { 255 /* log(u) = log( t u/t ) = log(t) + log(u/t) 256 log(t) is tabulated in the lookup table. 257 Express log(u/t) = log(1+z), where z = u/t - 1 = (u-t)/t. 258 cf. Cody & Waite. */ 259 z = (u.value - t.value) / t.value; 260 } 261 /* Series expansion of log(1+z). */ 262 w = z * z; 263 y = ((((((((((((l15 * z 264 + l14) * z 265 + l13) * z 266 + l12) * z 267 + l11) * z 268 + l10) * z 269 + l9) * z 270 + l8) * z 271 + l7) * z 272 + l6) * z 273 + l5) * z 274 + l4) * z 275 + l3) * z * w; 276 y -= 0.5 * w; 277 y += e * ln2b; /* Base 2 exponent offset times ln(2). */ 278 y += z; 279 y += logtbl[k-26]; /* log(t) - (t-1) */ 280 y += (t.value - 1.0L); 281 y += e * ln2a; 282 return y; 283 } 284