1141296Sdas 2141296Sdas/* @(#)e_log.c 1.3 95/01/18 */ 32116Sjkh/* 42116Sjkh * ==================================================== 52116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 62116Sjkh * 7141296Sdas * Developed at SunSoft, a Sun Microsystems, Inc. business. 82116Sjkh * Permission to use, copy, modify, and distribute this 9141296Sdas * software is freely granted, provided that this notice 102116Sjkh * is preserved. 112116Sjkh * ==================================================== 122116Sjkh */ 132116Sjkh 14176451Sdas#include <sys/cdefs.h> 15176451Sdas__FBSDID("$FreeBSD: releng/10.2/lib/msun/src/k_log.h 226376 2011-10-15 05:23:28Z das $"); 162116Sjkh 17226376Sdas/* 18226376Sdas * k_log1p(f): 19226376Sdas * Return log(1+f) - f for 1+f in ~[sqrt(2)/2, sqrt(2)]. 202116Sjkh * 21216210Sdas * The following describes the overall strategy for computing 22216210Sdas * logarithms in base e. The argument reduction and adding the final 23216210Sdas * term of the polynomial are done by the caller for increased accuracy 24216210Sdas * when different bases are used. 25216210Sdas * 26141296Sdas * Method : 27141296Sdas * 1. Argument Reduction: find k and f such that 28141296Sdas * x = 2^k * (1+f), 292116Sjkh * where sqrt(2)/2 < 1+f < sqrt(2) . 302116Sjkh * 312116Sjkh * 2. Approximation of log(1+f). 322116Sjkh * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) 332116Sjkh * = 2s + 2/3 s**3 + 2/5 s**5 + ....., 342116Sjkh * = 2s + s*R 35141296Sdas * We use a special Reme algorithm on [0,0.1716] to generate 36141296Sdas * a polynomial of degree 14 to approximate R The maximum error 372116Sjkh * of this polynomial approximation is bounded by 2**-58.45. In 382116Sjkh * other words, 392116Sjkh * 2 4 6 8 10 12 14 402116Sjkh * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s 412116Sjkh * (the values of Lg1 to Lg7 are listed in the program) 422116Sjkh * and 432116Sjkh * | 2 14 | -58.45 44141296Sdas * | Lg1*s +...+Lg7*s - R(z) | <= 2 452116Sjkh * | | 462116Sjkh * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. 472116Sjkh * In order to guarantee error in log below 1ulp, we compute log 482116Sjkh * by 492116Sjkh * log(1+f) = f - s*(f - R) (if f is not too large) 502116Sjkh * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) 51141296Sdas * 52141296Sdas * 3. Finally, log(x) = k*ln2 + log(1+f). 532116Sjkh * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) 54141296Sdas * Here ln2 is split into two floating point number: 552116Sjkh * ln2_hi + ln2_lo, 562116Sjkh * where n*ln2_hi is always exact for |n| < 2000. 572116Sjkh * 582116Sjkh * Special cases: 59141296Sdas * log(x) is NaN with signal if x < 0 (including -INF) ; 602116Sjkh * log(+INF) is +INF; log(0) is -INF with signal; 612116Sjkh * log(NaN) is that NaN with no signal. 622116Sjkh * 632116Sjkh * Accuracy: 642116Sjkh * according to an error analysis, the error is always less than 652116Sjkh * 1 ulp (unit in the last place). 662116Sjkh * 672116Sjkh * Constants: 68141296Sdas * The hexadecimal values are the intended ones for the following 69141296Sdas * constants. The decimal values may be used, provided that the 70141296Sdas * compiler will convert from decimal to binary accurately enough 712116Sjkh * to produce the hexadecimal values shown. 722116Sjkh */ 732116Sjkh 742116Sjkhstatic const double 752116SjkhLg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ 762116SjkhLg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ 772116SjkhLg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ 782116SjkhLg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ 792116SjkhLg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ 802116SjkhLg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ 812116SjkhLg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ 822116Sjkh 83216210Sdas/* 84226376Sdas * We always inline k_log1p(), since doing so produces a 85216210Sdas * substantial performance improvement (~40% on amd64). 86216210Sdas */ 87216210Sdasstatic inline double 88226376Sdask_log1p(double f) 892116Sjkh{ 90226376Sdas double hfsq,s,z,R,w,t1,t2; 912116Sjkh 92226376Sdas s = f/(2.0+f); 932116Sjkh z = s*s; 942116Sjkh w = z*z; 95226376Sdas t1= w*(Lg2+w*(Lg4+w*Lg6)); 96226376Sdas t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); 972116Sjkh R = t2+t1; 98226376Sdas hfsq=0.5*f*f; 99226376Sdas return s*(hfsq+R); 1002116Sjkh} 101