1135446Strhodes/* logll.c 2262706Serwin * 3135446Strhodes * Natural logarithm for 128-bit long double precision. 4135446Strhodes * 5174187Sdougb * 6135446Strhodes * 7135446Strhodes * SYNOPSIS: 8135446Strhodes * 9135446Strhodes * long double x, y, logq(); 10135446Strhodes * 11135446Strhodes * y = logq( x ); 12135446Strhodes * 13135446Strhodes * 14135446Strhodes * 15135446Strhodes * DESCRIPTION: 16135446Strhodes * 17135446Strhodes * Returns the base e (2.718...) logarithm of x. 18135446Strhodes * 19234010Sdougb * The argument is separated into its exponent and fractional 20135446Strhodes * parts. Use of a lookup table increases the speed of the routine. 21170222Sdougb * The program uses logarithms tabulated at intervals of 1/128 to 22135446Strhodes * cover the domain from approximately 0.7 to 1.4. 23135446Strhodes * 24135446Strhodes * On the interval [-1/128, +1/128] the logarithm of 1+x is approximated by 25135446Strhodes * log(1+x) = x - 0.5 x^2 + x^3 P(x) . 26135446Strhodes * 27135446Strhodes * 28135446Strhodes * 29218384Sdougb * ACCURACY: 30135446Strhodes * 31135446Strhodes * Relative error: 32193149Sdougb * arithmetic domain # trials peak rms 33135446Strhodes * IEEE 0.875, 1.125 100000 1.2e-34 4.1e-35 34135446Strhodes * IEEE 0.125, 8 100000 1.2e-34 4.1e-35 35135446Strhodes * 36135446Strhodes * 37193149Sdougb * WARNING: 38135446Strhodes * 39135446Strhodes * This program uses integer operations on bit fields of floating-point 40135446Strhodes * numbers. It does not work with data structures other than the 41135446Strhodes * structure assumed. 42135446Strhodes * 43135446Strhodes */ 44135446Strhodes 45135446Strhodes/* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov> 46135446Strhodes 47135446Strhodes This library is free software; you can redistribute it and/or 48135446Strhodes modify it under the terms of the GNU Lesser General Public 49135446Strhodes License as published by the Free Software Foundation; either 50135446Strhodes version 2.1 of the License, or (at your option) any later version. 51218384Sdougb 52218384Sdougb This library is distributed in the hope that it will be useful, 53218384Sdougb but WITHOUT ANY WARRANTY; without even the implied warranty of 54218384Sdougb MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 55135446Strhodes Lesser General Public License for more details. 56135446Strhodes 57135446Strhodes You should have received a copy of the GNU Lesser General Public 58170222Sdougb License along with this library; if not, see 59170222Sdougb <http://www.gnu.org/licenses/>. */ 60170222Sdougb 61170222Sdougb#include "quadmath-imp.h" 62170222Sdougb 63135446Strhodes/* log(1+x) = x - .5 x^2 + x^3 l(x) 64135446Strhodes -.0078125 <= x <= +.0078125 65135446Strhodes peak relative error 1.2e-37 */ 66135446Strhodesstatic const __float128 67135446Strhodesl3 = 3.333333333333333333333333333333336096926E-1Q, 68135446Strhodesl4 = -2.499999999999999999999999999486853077002E-1Q, 69135446Strhodesl5 = 1.999999999999999999999999998515277861905E-1Q, 70135446Strhodesl6 = -1.666666666666666666666798448356171665678E-1Q, 71135446Strhodesl7 = 1.428571428571428571428808945895490721564E-1Q, 72135446Strhodesl8 = -1.249999999999999987884655626377588149000E-1Q, 73135446Strhodesl9 = 1.111111111111111093947834982832456459186E-1Q, 74135446Strhodesl10 = -1.000000000000532974938900317952530453248E-1Q, 75135446Strhodesl11 = 9.090909090915566247008015301349979892689E-2Q, 76135446Strhodesl12 = -8.333333211818065121250921925397567745734E-2Q, 77135446Strhodesl13 = 7.692307559897661630807048686258659316091E-2Q, 78135446Strhodesl14 = -7.144242754190814657241902218399056829264E-2Q, 79135446Strhodesl15 = 6.668057591071739754844678883223432347481E-2Q; 80135446Strhodes 81135446Strhodes/* Lookup table of ln(t) - (t-1) 82135446Strhodes t = 0.5 + (k+26)/128) 83135446Strhodes k = 0, ..., 91 */ 84135446Strhodesstatic const __float128 logtbl[92] = { 85135446Strhodes-5.5345593589352099112142921677820359632418E-2Q, 86135446Strhodes-5.2108257402767124761784665198737642086148E-2Q, 87135446Strhodes-4.8991686870576856279407775480686721935120E-2Q, 88135446Strhodes-4.5993270766361228596215288742353061431071E-2Q, 89135446Strhodes-4.3110481649613269682442058976885699556950E-2Q, 90135446Strhodes-4.0340872319076331310838085093194799765520E-2Q, 91135446Strhodes-3.7682072451780927439219005993827431503510E-2Q, 92135446Strhodes-3.5131785416234343803903228503274262719586E-2Q, 93135446Strhodes-3.2687785249045246292687241862699949178831E-2Q, 94193149Sdougb-3.0347913785027239068190798397055267411813E-2Q, 95193149Sdougb-2.8110077931525797884641940838507561326298E-2Q, 96193149Sdougb-2.5972247078357715036426583294246819637618E-2Q, 97193149Sdougb-2.3932450635346084858612873953407168217307E-2Q, 98135446Strhodes-2.1988775689981395152022535153795155900240E-2Q, 99135446Strhodes-2.0139364778244501615441044267387667496733E-2Q, 100135446Strhodes-1.8382413762093794819267536615342902718324E-2Q, 101135446Strhodes-1.6716169807550022358923589720001638093023E-2Q, 102135446Strhodes-1.5138929457710992616226033183958974965355E-2Q, 103135446Strhodes-1.3649036795397472900424896523305726435029E-2Q, 104135446Strhodes-1.2244881690473465543308397998034325468152E-2Q, 105135446Strhodes-1.0924898127200937840689817557742469105693E-2Q, 106135446Strhodes-9.6875626072830301572839422532631079809328E-3Q, 107135446Strhodes-8.5313926245226231463436209313499745894157E-3Q, 108135446Strhodes-7.4549452072765973384933565912143044991706E-3Q, 109135446Strhodes-6.4568155251217050991200599386801665681310E-3Q, 110170222Sdougb-5.5356355563671005131126851708522185605193E-3Q, 111170222Sdougb-4.6900728132525199028885749289712348829878E-3Q, 112170222Sdougb-3.9188291218610470766469347968659624282519E-3Q, 113170222Sdougb-3.2206394539524058873423550293617843896540E-3Q, 114186462Sdougb-2.5942708080877805657374888909297113032132E-3Q, 115186462Sdougb-2.0385211375711716729239156839929281289086E-3Q, 116186462Sdougb-1.5522183228760777967376942769773768850872E-3Q, 117186462Sdougb-1.1342191863606077520036253234446621373191E-3Q, 118186462Sdougb-7.8340854719967065861624024730268350459991E-4Q, 119186462Sdougb-4.9869831458030115699628274852562992756174E-4Q, 120170222Sdougb-2.7902661731604211834685052867305795169688E-4Q, 121170222Sdougb-1.2335696813916860754951146082826952093496E-4Q, 122170222Sdougb-3.0677461025892873184042490943581654591817E-5Q, 123170222Sdougb#define ZERO logtbl[38] 124170222Sdougb 0.0000000000000000000000000000000000000000E0Q, 125170222Sdougb-3.0359557945051052537099938863236321874198E-5Q, 126170222Sdougb-1.2081346403474584914595395755316412213151E-4Q, 127170222Sdougb-2.7044071846562177120083903771008342059094E-4Q, 128186462Sdougb-4.7834133324631162897179240322783590830326E-4Q, 129186462Sdougb-7.4363569786340080624467487620270965403695E-4Q, 130186462Sdougb-1.0654639687057968333207323853366578860679E-3Q, 131186462Sdougb-1.4429854811877171341298062134712230604279E-3Q, 132186462Sdougb-1.8753781835651574193938679595797367137975E-3Q, 133186462Sdougb-2.3618380914922506054347222273705859653658E-3Q, 134170222Sdougb-2.9015787624124743013946600163375853631299E-3Q, 135170222Sdougb-3.4938307889254087318399313316921940859043E-3Q, 136170222Sdougb-4.1378413103128673800485306215154712148146E-3Q, 137170222Sdougb-4.8328735414488877044289435125365629849599E-3Q, 138170222Sdougb-5.5782063183564351739381962360253116934243E-3Q, 139170222Sdougb-6.3731336597098858051938306767880719015261E-3Q, 140170222Sdougb-7.2169643436165454612058905294782949315193E-3Q, 141170222Sdougb-8.1090214990427641365934846191367315083867E-3Q, 142186462Sdougb-9.0486422112807274112838713105168375482480E-3Q, 143186462Sdougb-1.0035177140880864314674126398350812606841E-2Q, 144186462Sdougb-1.1067990155502102718064936259435676477423E-2Q, 145186462Sdougb-1.2146457974158024928196575103115488672416E-2Q, 146186462Sdougb-1.3269969823361415906628825374158424754308E-2Q, 147186462Sdougb-1.4437927104692837124388550722759686270765E-2Q, 148170222Sdougb-1.5649743073340777659901053944852735064621E-2Q, 149170222Sdougb-1.6904842527181702880599758489058031645317E-2Q, 150170222Sdougb-1.8202661505988007336096407340750378994209E-2Q, 151170222Sdougb-1.9542647000370545390701192438691126552961E-2Q, 152170222Sdougb-2.0924256670080119637427928803038530924742E-2Q, 153170222Sdougb-2.2346958571309108496179613803760727786257E-2Q, 154170222Sdougb-2.3810230892650362330447187267648486279460E-2Q, 155170222Sdougb-2.5313561699385640380910474255652501521033E-2Q, 156186462Sdougb-2.6856448685790244233704909690165496625399E-2Q, 157186462Sdougb-2.8438398935154170008519274953860128449036E-2Q, 158186462Sdougb-3.0058928687233090922411781058956589863039E-2Q, 159186462Sdougb-3.1717563112854831855692484086486099896614E-2Q, 160186462Sdougb-3.3413836095418743219397234253475252001090E-2Q, 161186462Sdougb-3.5147290019036555862676702093393332533702E-2Q, 162170222Sdougb-3.6917475563073933027920505457688955423688E-2Q, 163170222Sdougb-3.8723951502862058660874073462456610731178E-2Q, 164170222Sdougb-4.0566284516358241168330505467000838017425E-2Q, 165170222Sdougb-4.2444048996543693813649967076598766917965E-2Q, 166170222Sdougb-4.4356826869355401653098777649745233339196E-2Q, 167170222Sdougb-4.6304207416957323121106944474331029996141E-2Q, 168170222Sdougb-4.8285787106164123613318093945035804818364E-2Q, 169170222Sdougb-5.0301169421838218987124461766244507342648E-2Q, 170186462Sdougb-5.2349964705088137924875459464622098310997E-2Q, 171186462Sdougb-5.4431789996103111613753440311680967840214E-2Q, 172186462Sdougb-5.6546268881465384189752786409400404404794E-2Q, 173186462Sdougb-5.8693031345788023909329239565012647817664E-2Q, 174186462Sdougb-6.0871713627532018185577188079210189048340E-2Q, 175186462Sdougb-6.3081958078862169742820420185833800925568E-2Q, 176170222Sdougb-6.5323413029406789694910800219643791556918E-2Q, 177170222Sdougb-6.7595732653791419081537811574227049288168E-2Q 178170222Sdougb}; 179170222Sdougb 180135446Strhodes/* ln(2) = ln2a + ln2b with extended precision. */ 181135446Strhodesstatic const __float128 182135446Strhodes ln2a = 6.93145751953125e-1Q, 183135446Strhodes ln2b = 1.4286068203094172321214581765680755001344E-6Q; 184135446Strhodes 185135446Strhodes__float128 186135446Strhodeslogq(__float128 x) 187135446Strhodes{ 188193149Sdougb __float128 z, y, w; 189193149Sdougb ieee854_float128 u, t; 190193149Sdougb unsigned int m; 191193149Sdougb int k, e; 192193149Sdougb 193135446Strhodes u.value = x; 194135446Strhodes m = u.words32.w0; 195135446Strhodes 196135446Strhodes /* Check for IEEE special cases. */ 197193149Sdougb k = m & 0x7fffffff; 198135446Strhodes /* log(0) = -infinity. */ 199135446Strhodes if ((k | u.words32.w1 | u.words32.w2 | u.words32.w3) == 0) 200135446Strhodes { 201135446Strhodes return -0.5Q / ZERO; 202135446Strhodes } 203135446Strhodes /* log ( x < 0 ) = NaN */ 204135446Strhodes if (m & 0x80000000) 205193149Sdougb { 206224092Sdougb return (x - x) / ZERO; 207193149Sdougb } 208224092Sdougb /* log (infinity or NaN) */ 209224092Sdougb if (k >= 0x7fff0000) 210193149Sdougb { 211135446Strhodes return x + x; 212135446Strhodes } 213135446Strhodes 214193149Sdougb /* Extract exponent and reduce domain to 0.703125 <= u < 1.40625 */ 215193149Sdougb u.value = frexpq (x, &e); 216193149Sdougb m = u.words32.w0 & 0xffff; 217193149Sdougb m |= 0x10000; 218193149Sdougb /* Find lookup table index k from high order bits of the significand. */ 219193149Sdougb if (m < 0x16800) 220193149Sdougb { 221193149Sdougb k = (m - 0xff00) >> 9; 222193149Sdougb /* t is the argument 0.5 + (k+26)/128 223135446Strhodes of the nearest item to u in the lookup table. */ 224135446Strhodes t.words32.w0 = 0x3fff0000 + (k << 9); 225224092Sdougb t.words32.w1 = 0; 226224092Sdougb t.words32.w2 = 0; 227224092Sdougb t.words32.w3 = 0; 228224092Sdougb u.words32.w0 += 0x10000; 229224092Sdougb e -= 1; 230224092Sdougb k += 64; 231224092Sdougb } 232224092Sdougb else 233224092Sdougb { 234224092Sdougb k = (m - 0xfe00) >> 10; 235224092Sdougb t.words32.w0 = 0x3ffe0000 + (k << 10); 236224092Sdougb t.words32.w1 = 0; 237224092Sdougb t.words32.w2 = 0; 238224092Sdougb t.words32.w3 = 0; 239224092Sdougb } 240224092Sdougb /* On this interval the table is not used due to cancellation error. */ 241224092Sdougb if ((x <= 1.0078125Q) && (x >= 0.9921875Q)) 242254402Serwin { 243254402Serwin if (x == 1) 244254402Serwin return 0; 245224092Sdougb z = x - 1; 246224092Sdougb k = 64; 247224092Sdougb t.value = 1; 248224092Sdougb e = 0; 249224092Sdougb } 250224092Sdougb else 251224092Sdougb { 252224092Sdougb /* log(u) = log( t u/t ) = log(t) + log(u/t) 253224092Sdougb log(t) is tabulated in the lookup table. 254224092Sdougb Express log(u/t) = log(1+z), where z = u/t - 1 = (u-t)/t. 255224092Sdougb cf. Cody & Waite. */ 256224092Sdougb z = (u.value - t.value) / t.value; 257224092Sdougb } 258224092Sdougb /* Series expansion of log(1+z). */ 259224092Sdougb w = z * z; 260224092Sdougb y = ((((((((((((l15 * z 261224092Sdougb + l14) * z 262224092Sdougb + l13) * z 263224092Sdougb + l12) * z 264224092Sdougb + l11) * z 265224092Sdougb + l10) * z 266224092Sdougb + l9) * z 267224092Sdougb + l8) * z 268224092Sdougb + l7) * z 269224092Sdougb + l6) * z 270224092Sdougb + l5) * z 271224092Sdougb + l4) * z 272224092Sdougb + l3) * z * w; 273224092Sdougb y -= 0.5 * w; 274224092Sdougb y += e * ln2b; /* Base 2 exponent offset times ln(2). */ 275224092Sdougb y += z; 276224092Sdougb y += logtbl[k-26]; /* log(t) - (t-1) */ 277224092Sdougb y += (t.value - 1); 278224092Sdougb y += e * ln2a; 279224092Sdougb return y; 280224092Sdougb} 281224092Sdougb