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 */ 66static const __float128 67l3 = 3.333333333333333333333333333333336096926E-1Q, 68l4 = -2.499999999999999999999999999486853077002E-1Q, 69l5 = 1.999999999999999999999999998515277861905E-1Q, 70l6 = -1.666666666666666666666798448356171665678E-1Q, 71l7 = 1.428571428571428571428808945895490721564E-1Q, 72l8 = -1.249999999999999987884655626377588149000E-1Q, 73l9 = 1.111111111111111093947834982832456459186E-1Q, 74l10 = -1.000000000000532974938900317952530453248E-1Q, 75l11 = 9.090909090915566247008015301349979892689E-2Q, 76l12 = -8.333333211818065121250921925397567745734E-2Q, 77l13 = 7.692307559897661630807048686258659316091E-2Q, 78l14 = -7.144242754190814657241902218399056829264E-2Q, 79l15 = 6.668057591071739754844678883223432347481E-2Q; 80 81/* Lookup table of ln(t) - (t-1) 82 t = 0.5 + (k+26)/128) 83 k = 0, ..., 91 */ 84static 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. */ 181static const __float128 182 ln2a = 6.93145751953125e-1Q, 183 ln2b = 1.4286068203094172321214581765680755001344E-6Q; 184 185__float128 186logq(__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} 281