1/* $NetBSD$ */ 2 3/* 4 * Copyright (c) 1992, 1993 5 * The Regents of the University of California. All rights reserved. 6 * 7 * Redistribution and use in source and binary forms, with or without 8 * modification, are permitted provided that the following conditions 9 * are met: 10 * 1. Redistributions of source code must retain the above copyright 11 * notice, this list of conditions and the following disclaimer. 12 * 2. Redistributions in binary form must reproduce the above copyright 13 * notice, this list of conditions and the following disclaimer in the 14 * documentation and/or other materials provided with the distribution. 15 * 3. All advertising materials mentioning features or use of this software 16 * must display the following acknowledgement: 17 * This product includes software developed by the University of 18 * California, Berkeley and its contributors. 19 * 4. Neither the name of the University nor the names of its contributors 20 * may be used to endorse or promote products derived from this software 21 * without specific prior written permission. 22 * 23 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 24 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 25 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 26 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 27 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 28 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 29 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 30 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 31 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 32 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 33 * SUCH DAMAGE. 34 */ 35 36/* @(#)log.c 8.2 (Berkeley) 11/30/93 */ 37#include <sys/cdefs.h> 38#if 0 39__FBSDID("$FreeBSD: release/9.0.0/lib/msun/bsdsrc/b_log.c 176449 2008-02-22 02:26:51Z das $"); 40#else 41__RCSID("$NetBSD$"); 42#endif 43 44#include <errno.h> 45#include "math.h" 46#include "math_private.h" 47 48/* Table-driven natural logarithm. 49 * 50 * This code was derived, with minor modifications, from: 51 * Peter Tang, "Table-Driven Implementation of the 52 * Logarithm in IEEE Floating-Point arithmetic." ACM Trans. 53 * Math Software, vol 16. no 4, pp 378-400, Dec 1990). 54 * 55 * Calculates log(2^m*F*(1+f/F)), |f/j| <= 1/256, 56 * where F = j/128 for j an integer in [0, 128]. 57 * 58 * log(2^m) = log2_hi*m + log2_tail*m 59 * since m is an integer, the dominant term is exact. 60 * m has at most 10 digits (for subnormal numbers), 61 * and log2_hi has 11 trailing zero bits. 62 * 63 * log(F) = logF_hi[j] + logF_lo[j] is in tabular form in log_table.h 64 * logF_hi[] + 512 is exact. 65 * 66 * log(1+f/F) = 2*f/(2*F + f) + 1/12 * (2*f/(2*F + f))**3 + ... 67 * the leading term is calculated to extra precision in two 68 * parts, the larger of which adds exactly to the dominant 69 * m and F terms. 70 * There are two cases: 71 * 1. when m, j are non-zero (m | j), use absolute 72 * precision for the leading term. 73 * 2. when m = j = 0, |1-x| < 1/256, and log(x) ~= (x-1). 74 * In this case, use a relative precision of 24 bits. 75 * (This is done differently in the original paper) 76 * 77 * Special cases: 78 * 0 return signalling -Inf 79 * neg return signalling NaN 80 * +Inf return +Inf 81*/ 82 83#define N 128 84 85/* Table of log(Fj) = logF_head[j] + logF_tail[j], for Fj = 1+j/128. 86 * Used for generation of extend precision logarithms. 87 * The constant 35184372088832 is 2^45, so the divide is exact. 88 * It ensures correct reading of logF_head, even for inaccurate 89 * decimal-to-binary conversion routines. (Everybody gets the 90 * right answer for integers less than 2^53.) 91 * Values for log(F) were generated using error < 10^-57 absolute 92 * with the bc -l package. 93*/ 94static const double A1 = .08333333333333178827; 95static const double A2 = .01250000000377174923; 96static const double A3 = .002232139987919447809; 97static const double A4 = .0004348877777076145742; 98 99static const double logF_head[N+1] = { 100 0., 101 .007782140442060381246, 102 .015504186535963526694, 103 .023167059281547608406, 104 .030771658666765233647, 105 .038318864302141264488, 106 .045809536031242714670, 107 .053244514518837604555, 108 .060624621816486978786, 109 .067950661908525944454, 110 .075223421237524235039, 111 .082443669210988446138, 112 .089612158689760690322, 113 .096729626458454731618, 114 .103796793681567578460, 115 .110814366340264314203, 116 .117783035656430001836, 117 .124703478501032805070, 118 .131576357788617315236, 119 .138402322859292326029, 120 .145182009844575077295, 121 .151916042025732167530, 122 .158605030176659056451, 123 .165249572895390883786, 124 .171850256926518341060, 125 .178407657472689606947, 126 .184922338493834104156, 127 .191394852999565046047, 128 .197825743329758552135, 129 .204215541428766300668, 130 .210564769107350002741, 131 .216873938300523150246, 132 .223143551314024080056, 133 .229374101064877322642, 134 .235566071312860003672, 135 .241719936886966024758, 136 .247836163904594286577, 137 .253915209980732470285, 138 .259957524436686071567, 139 .265963548496984003577, 140 .271933715484010463114, 141 .277868451003087102435, 142 .283768173130738432519, 143 .289633292582948342896, 144 .295464212893421063199, 145 .301261330578199704177, 146 .307025035294827830512, 147 .312755710004239517729, 148 .318453731118097493890, 149 .324119468654316733591, 150 .329753286372579168528, 151 .335355541920762334484, 152 .340926586970454081892, 153 .346466767346100823488, 154 .351976423156884266063, 155 .357455888922231679316, 156 .362905493689140712376, 157 .368325561158599157352, 158 .373716409793814818840, 159 .379078352934811846353, 160 .384411698910298582632, 161 .389716751140440464951, 162 .394993808240542421117, 163 .400243164127459749579, 164 .405465108107819105498, 165 .410659924985338875558, 166 .415827895143593195825, 167 .420969294644237379543, 168 .426084395310681429691, 169 .431173464818130014464, 170 .436236766774527495726, 171 .441274560805140936281, 172 .446287102628048160113, 173 .451274644139630254358, 174 .456237433481874177232, 175 .461175715122408291790, 176 .466089729924533457960, 177 .470979715219073113985, 178 .475845904869856894947, 179 .480688529345570714212, 180 .485507815781602403149, 181 .490303988045525329653, 182 .495077266798034543171, 183 .499827869556611403822, 184 .504556010751912253908, 185 .509261901790523552335, 186 .513945751101346104405, 187 .518607764208354637958, 188 .523248143765158602036, 189 .527867089620485785417, 190 .532464798869114019908, 191 .537041465897345915436, 192 .541597282432121573947, 193 .546132437597407260909, 194 .550647117952394182793, 195 .555141507540611200965, 196 .559615787935399566777, 197 .564070138285387656651, 198 .568504735352689749561, 199 .572919753562018740922, 200 .577315365035246941260, 201 .581691739635061821900, 202 .586049045003164792433, 203 .590387446602107957005, 204 .594707107746216934174, 205 .599008189645246602594, 206 .603290851438941899687, 207 .607555250224322662688, 208 .611801541106615331955, 209 .616029877215623855590, 210 .620240409751204424537, 211 .624433288012369303032, 212 .628608659422752680256, 213 .632766669570628437213, 214 .636907462236194987781, 215 .641031179420679109171, 216 .645137961373620782978, 217 .649227946625615004450, 218 .653301272011958644725, 219 .657358072709030238911, 220 .661398482245203922502, 221 .665422632544505177065, 222 .669430653942981734871, 223 .673422675212350441142, 224 .677398823590920073911, 225 .681359224807238206267, 226 .685304003098281100392, 227 .689233281238557538017, 228 .693147180560117703862 229}; 230 231static const double logF_tail[N+1] = { 232 0., 233 -.00000000000000543229938420049, 234 .00000000000000172745674997061, 235 -.00000000000001323017818229233, 236 -.00000000000001154527628289872, 237 -.00000000000000466529469958300, 238 .00000000000005148849572685810, 239 -.00000000000002532168943117445, 240 -.00000000000005213620639136504, 241 -.00000000000001819506003016881, 242 .00000000000006329065958724544, 243 .00000000000008614512936087814, 244 -.00000000000007355770219435028, 245 .00000000000009638067658552277, 246 .00000000000007598636597194141, 247 .00000000000002579999128306990, 248 -.00000000000004654729747598444, 249 -.00000000000007556920687451336, 250 .00000000000010195735223708472, 251 -.00000000000017319034406422306, 252 -.00000000000007718001336828098, 253 .00000000000010980754099855238, 254 -.00000000000002047235780046195, 255 -.00000000000008372091099235912, 256 .00000000000014088127937111135, 257 .00000000000012869017157588257, 258 .00000000000017788850778198106, 259 .00000000000006440856150696891, 260 .00000000000016132822667240822, 261 -.00000000000007540916511956188, 262 -.00000000000000036507188831790, 263 .00000000000009120937249914984, 264 .00000000000018567570959796010, 265 -.00000000000003149265065191483, 266 -.00000000000009309459495196889, 267 .00000000000017914338601329117, 268 -.00000000000001302979717330866, 269 .00000000000023097385217586939, 270 .00000000000023999540484211737, 271 .00000000000015393776174455408, 272 -.00000000000036870428315837678, 273 .00000000000036920375082080089, 274 -.00000000000009383417223663699, 275 .00000000000009433398189512690, 276 .00000000000041481318704258568, 277 -.00000000000003792316480209314, 278 .00000000000008403156304792424, 279 -.00000000000034262934348285429, 280 .00000000000043712191957429145, 281 -.00000000000010475750058776541, 282 -.00000000000011118671389559323, 283 .00000000000037549577257259853, 284 .00000000000013912841212197565, 285 .00000000000010775743037572640, 286 .00000000000029391859187648000, 287 -.00000000000042790509060060774, 288 .00000000000022774076114039555, 289 .00000000000010849569622967912, 290 -.00000000000023073801945705758, 291 .00000000000015761203773969435, 292 .00000000000003345710269544082, 293 -.00000000000041525158063436123, 294 .00000000000032655698896907146, 295 -.00000000000044704265010452446, 296 .00000000000034527647952039772, 297 -.00000000000007048962392109746, 298 .00000000000011776978751369214, 299 -.00000000000010774341461609578, 300 .00000000000021863343293215910, 301 .00000000000024132639491333131, 302 .00000000000039057462209830700, 303 -.00000000000026570679203560751, 304 .00000000000037135141919592021, 305 -.00000000000017166921336082431, 306 -.00000000000028658285157914353, 307 -.00000000000023812542263446809, 308 .00000000000006576659768580062, 309 -.00000000000028210143846181267, 310 .00000000000010701931762114254, 311 .00000000000018119346366441110, 312 .00000000000009840465278232627, 313 -.00000000000033149150282752542, 314 -.00000000000018302857356041668, 315 -.00000000000016207400156744949, 316 .00000000000048303314949553201, 317 -.00000000000071560553172382115, 318 .00000000000088821239518571855, 319 -.00000000000030900580513238244, 320 -.00000000000061076551972851496, 321 .00000000000035659969663347830, 322 .00000000000035782396591276383, 323 -.00000000000046226087001544578, 324 .00000000000062279762917225156, 325 .00000000000072838947272065741, 326 .00000000000026809646615211673, 327 -.00000000000010960825046059278, 328 .00000000000002311949383800537, 329 -.00000000000058469058005299247, 330 -.00000000000002103748251144494, 331 -.00000000000023323182945587408, 332 -.00000000000042333694288141916, 333 -.00000000000043933937969737844, 334 .00000000000041341647073835565, 335 .00000000000006841763641591466, 336 .00000000000047585534004430641, 337 .00000000000083679678674757695, 338 -.00000000000085763734646658640, 339 .00000000000021913281229340092, 340 -.00000000000062242842536431148, 341 -.00000000000010983594325438430, 342 .00000000000065310431377633651, 343 -.00000000000047580199021710769, 344 -.00000000000037854251265457040, 345 .00000000000040939233218678664, 346 .00000000000087424383914858291, 347 .00000000000025218188456842882, 348 -.00000000000003608131360422557, 349 -.00000000000050518555924280902, 350 .00000000000078699403323355317, 351 -.00000000000067020876961949060, 352 .00000000000016108575753932458, 353 .00000000000058527188436251509, 354 -.00000000000035246757297904791, 355 -.00000000000018372084495629058, 356 .00000000000088606689813494916, 357 .00000000000066486268071468700, 358 .00000000000063831615170646519, 359 .00000000000025144230728376072, 360 -.00000000000017239444525614834 361}; 362 363/* 364 * Extra precision variant, returning struct {double a, b;}; 365 * log(x) = a+b to 63 bits, with a rounded to 26 bits. 366 */ 367struct Double 368__log__D(double x) 369{ 370 int m, j; 371 double F, f, g, q, u, v, u2; 372 volatile double u1; 373 struct Double r; 374 375 /* Argument reduction: 1 <= g < 2; x/2^m = g; */ 376 /* y = F*(1 + f/F) for |f| <= 2^-8 */ 377 378 m = logb(x); 379 g = ldexp(x, -m); 380 if (m == -1022) { 381 j = logb(g), m += j; 382 g = ldexp(g, -j); 383 } 384 j = N*(g-1) + .5; 385 F = (1.0/N) * j + 1; 386 f = g - F; 387 388 g = 1/(2*F+f); 389 u = 2*f*g; 390 v = u*u; 391 q = u*v*(A1 + v*(A2 + v*(A3 + v*A4))); 392 if (m | j) 393 u1 = u + 513, u1 -= 513; 394 else 395 u1 = u, TRUNC(u1); 396 u2 = (2.0*(f - F*u1) - u1*f) * g; 397 398 u1 += m*logF_head[N] + logF_head[j]; 399 400 u2 += logF_tail[j]; u2 += q; 401 u2 += logF_tail[N]*m; 402 r.a = u1 + u2; /* Only difference is here */ 403 TRUNC(r.a); 404 r.b = (u1 - r.a) + u2; 405 return (r); 406} 407