diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/e_log.c')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_log.c | 348 |
1 files changed, 190 insertions, 158 deletions
diff --git a/sysdeps/ieee754/dbl-64/e_log.c b/sysdeps/ieee754/dbl-64/e_log.c index 851bd30198..e55d74e561 100644 --- a/sysdeps/ieee754/dbl-64/e_log.c +++ b/sysdeps/ieee754/dbl-64/e_log.c @@ -1,165 +1,197 @@ -/* @(#)e_log.c 5.1 93/09/24 */ /* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * IBM Accurate Mathematical Library + * Copyright (c) International Business Machines Corp., 2001 * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ -/* Modified by Naohiko Shimizu/Tokai University, Japan 1997/08/25, - for performance improvement on pipelined processors. -*/ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_log.c,v 1.8 1995/05/10 20:45:49 jtc Exp $"; -#endif - -/* __ieee754_log(x) - * Return the logarithm of x - * - * Method : - * 1. Argument Reduction: find k and f such that - * x = 2^k * (1+f), - * where sqrt(2)/2 < 1+f < sqrt(2) . - * - * 2. Approximation of log(1+f). - * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) - * = 2s + 2/3 s**3 + 2/5 s**5 + ....., - * = 2s + s*R - * We use a special Reme algorithm on [0,0.1716] to generate - * a polynomial of degree 14 to approximate R The maximum error - * of this polynomial approximation is bounded by 2**-58.45. In - * other words, - * 2 4 6 8 10 12 14 - * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s - * (the values of Lg1 to Lg7 are listed in the program) - * and - * | 2 14 | -58.45 - * | Lg1*s +...+Lg7*s - R(z) | <= 2 - * | | - * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. - * In order to guarantee error in log below 1ulp, we compute log - * by - * log(1+f) = f - s*(f - R) (if f is not too large) - * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. * - * 3. Finally, log(x) = k*ln2 + log(1+f). - * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) - * Here ln2 is split into two floating point number: - * ln2_hi + ln2_lo, - * where n*ln2_hi is always exact for |n| < 2000. + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. * - * Special cases: - * log(x) is NaN with signal if x < 0 (including -INF) ; - * log(+INF) is +INF; log(0) is -INF with signal; - * log(NaN) is that NaN with no signal. - * - * Accuracy: - * according to an error analysis, the error is always less than - * 1 ulp (unit in the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the - * compiler will convert from decimal to binary accurately enough - * to produce the hexadecimal values shown. + * You should have received a copy of the GNU Lesser General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ +/*********************************************************************/ +/* */ +/* MODULE_NAME:ulog.h */ +/* */ +/* FUNCTION:ulog */ +/* */ +/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h ulog.h */ +/* mpexp.c mplog.c mpa.c */ +/* ulog.tbl */ +/* */ +/* An ultimate log routine. Given an IEEE double machine number x */ +/* it computes the correctly rounded (to nearest) value of log(x). */ +/* Assumption: Machine arithmetic operations are performed in */ +/* round to nearest mode of IEEE 754 standard. */ +/* */ +/*********************************************************************/ + + +#include "endian.h" +#include "dla.h" +#include "mpa.h" +#include "MathLib.h" +void __mplog(mp_no *, mp_no *, int); + +/*********************************************************************/ +/* An ultimate log routine. Given an IEEE double machine number x */ +/* it computes the correctly rounded (to nearest) value of log(x). */ +/*********************************************************************/ +double __ieee754_log(double x) { +#define M 4 + static const int pr[M]={8,10,18,32}; + int i,j,k,n,ux,dx,p; + double dbl_n,u,p0,q,r0,w,nln2a,luai,lubi,lvaj,lvbj, + sij,ssij,ttij,A,B,B0,y,y1,y2,polI,polII,sa,sb, + t1,t2,t3,t4,t5,t6,t7,t8,t,ra,rb,ww, + a0,aa0,s1,s2,ss2,s3,ss3,a1,aa1,a,aa,b,bb,c; + number num; + mp_no mpx,mpy,mpy1,mpy2,mperr; + +#include "ulog.tbl" +#include "ulog.h" + + /* Treating special values of x ( x<=0, x=INF, x=NaN etc.). */ + + num.d = x; ux = num.i[HIGH_HALF]; dx = num.i[LOW_HALF]; + n=0; + if (ux < 0x00100000) { + if (((ux & 0x7fffffff) | dx) == 0) return MHALF/ZERO; /* return -INF */ + if (ux < 0) return (x-x)/ZERO; /* return NaN */ + n -= 54; x *= two54.d; /* scale x */ + num.d = x; + } + if (ux >= 0x7ff00000) return x+x; /* INF or NaN */ + + /* Regular values of x */ + + w = x-ONE; + if (ABS(w) > U03) { goto case_03; } + + + /*--- Stage I, the case abs(x-1) < 0.03 */ + + t8 = MHALF*w; + EMULV(t8,w,a,aa,t1,t2,t3,t4,t5) + EADD(w,a,b,bb) + + /* Evaluate polynomial II */ + polII = (b0.d+w*(b1.d+w*(b2.d+w*(b3.d+w*(b4.d+ + w*(b5.d+w*(b6.d+w*(b7.d+w*b8.d))))))))*w*w*w; + c = (aa+bb)+polII; + + /* End stage I, case abs(x-1) < 0.03 */ + if ((y=b+(c+b*E2)) == b+(c-b*E2)) return y; + + /*--- Stage II, the case abs(x-1) < 0.03 */ + + a = d11.d+w*(d12.d+w*(d13.d+w*(d14.d+w*(d15.d+w*(d16.d+ + w*(d17.d+w*(d18.d+w*(d19.d+w*d20.d)))))))); + EMULV(w,a,s2,ss2,t1,t2,t3,t4,t5) + ADD2(d10.d,dd10.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d9.d,dd9.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d8.d,dd8.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d7.d,dd7.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d6.d,dd6.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d5.d,dd5.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d4.d,dd4.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d3.d,dd3.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d2.d,dd2.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + MUL2(w,ZERO,s2,ss2,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(w,ZERO, s3,ss3, b, bb,t1,t2) + + /* End stage II, case abs(x-1) < 0.03 */ + if ((y=b+(bb+b*E4)) == b+(bb-b*E4)) return y; + goto stage_n; + + /*--- Stage I, the case abs(x-1) > 0.03 */ + case_03: + + /* Find n,u such that x = u*2**n, 1/sqrt(2) < u < sqrt(2) */ + n += (num.i[HIGH_HALF] >> 20) - 1023; + num.i[HIGH_HALF] = (num.i[HIGH_HALF] & 0x000fffff) | 0x3ff00000; + if (num.d > SQRT_2) { num.d *= HALF; n++; } + u = num.d; dbl_n = (double) n; + + /* Find i such that ui=1+(i-75)/2**8 is closest to u (i= 0,1,2,...,181) */ + num.d += h1.d; + i = (num.i[HIGH_HALF] & 0x000fffff) >> 12; + + /* Find j such that vj=1+(j-180)/2**16 is closest to v=u/ui (j= 0,...,361) */ + num.d = u*Iu[i].d + h2.d; + j = (num.i[HIGH_HALF] & 0x000fffff) >> 4; + + /* Compute w=(u-ui*vj)/(ui*vj) */ + p0=(ONE+(i-75)*DEL_U)*(ONE+(j-180)*DEL_V); + q=u-p0; r0=Iu[i].d*Iv[j].d; w=q*r0; + + /* Evaluate polynomial I */ + polI = w+(a2.d+a3.d*w)*w*w; + + /* Add up everything */ + nln2a = dbl_n*LN2A; + luai = Lu[i][0].d; lubi = Lu[i][1].d; + lvaj = Lv[j][0].d; lvbj = Lv[j][1].d; + EADD(luai,lvaj,sij,ssij) + EADD(nln2a,sij,A ,ttij) + B0 = (((lubi+lvbj)+ssij)+ttij)+dbl_n*LN2B; + B = polI+B0; + + /* End stage I, case abs(x-1) >= 0.03 */ + if ((y=A+(B+E1)) == A+(B-E1)) return y; + + + /*--- Stage II, the case abs(x-1) > 0.03 */ + + /* Improve the accuracy of r0 */ + EMULV(p0,r0,sa,sb,t1,t2,t3,t4,t5) + t=r0*((ONE-sa)-sb); + EADD(r0,t,ra,rb) + + /* Compute w */ + MUL2(q,ZERO,ra,rb,w,ww,t1,t2,t3,t4,t5,t6,t7,t8) + + EADD(A,B0,a0,aa0) + + /* Evaluate polynomial III */ + s1 = (c3.d+(c4.d+c5.d*w)*w)*w; + EADD(c2.d,s1,s2,ss2) + MUL2(s2,ss2,w,ww,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8) + MUL2(s3,ss3,w,ww,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(s2,ss2,w,ww,s3,ss3,t1,t2) + ADD2(s3,ss3,a0,aa0,a1,aa1,t1,t2) + + /* End stage II, case abs(x-1) >= 0.03 */ + if ((y=a1+(aa1+E3)) == a1+(aa1-E3)) return y; + + + /* Final stages. Use multi-precision arithmetic. */ + stage_n: -#include "math.h" -#include "math_private.h" -#define half Lg[8] -#define two Lg[9] -#ifdef __STDC__ -static const double -#else -static double -#endif -ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ -ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ -two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ - Lg[] = {0.0, - 6.666666666666735130e-01, /* 3FE55555 55555593 */ - 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ - 2.857142874366239149e-01, /* 3FD24924 94229359 */ - 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ - 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ - 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ - 1.479819860511658591e-01, /* 3FC2F112 DF3E5244 */ - 0.5, - 2.0}; -#ifdef __STDC__ -static const double zero = 0.0; -#else -static double zero = 0.0; -#endif - -#ifdef __STDC__ - double __ieee754_log(double x) -#else - double __ieee754_log(x) - double x; -#endif -{ - double hfsq,f,s,z,R,w,dk,t11,t12,t21,t22,w2,zw2; -#ifdef DO_NOT_USE_THIS - double t1,t2; -#endif - int32_t k,hx,i,j; - u_int32_t lx; - - EXTRACT_WORDS(hx,lx,x); - - k=0; - if (hx < 0x00100000) { /* x < 2**-1022 */ - if (((hx&0x7fffffff)|lx)==0) - return -two54/(x-x); /* log(+-0)=-inf */ - if (hx<0) return (x-x)/(x-x); /* log(-#) = NaN */ - k -= 54; x *= two54; /* subnormal number, scale up x */ - GET_HIGH_WORD(hx,x); - } - if (hx >= 0x7ff00000) return x+x; - k += (hx>>20)-1023; - hx &= 0x000fffff; - i = (hx+0x95f64)&0x100000; - SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */ - k += (i>>20); - f = x-1.0; - if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */ - if(f==zero) { - if(k==0) return zero; else {dk=(double)k; - return dk*ln2_hi+dk*ln2_lo;} - } - R = f*f*(half-0.33333333333333333*f); - if(k==0) return f-R; else {dk=(double)k; - return dk*ln2_hi-((R-dk*ln2_lo)-f);} - } - s = f/(two+f); - dk = (double)k; - z = s*s; - i = hx-0x6147a; - w = z*z; - j = 0x6b851-hx; -#ifdef DO_NOT_USE_THIS - t1= w*(Lg2+w*(Lg4+w*Lg6)); - t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); - R = t2+t1; -#else - t21 = Lg[5]+w*Lg[7]; w2=w*w; - t22 = Lg[1]+w*Lg[3]; zw2=z*w2; - t11 = Lg[4]+w*Lg[6]; - t12 = w*Lg[2]; - R = t12 + w2*t11 + z*t22 + zw2*t21; -#endif - i |= j; - if(i>0) { - hfsq=0.5*f*f; - if(k==0) return f-(hfsq-s*(hfsq+R)); else - return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f); - } else { - if(k==0) return f-s*(f-R); else - return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f); - } + for (i=0; i<M; i++) { + p = pr[i]; + dbl_mp(x,&mpx,p); dbl_mp(y,&mpy,p); + __mplog(&mpx,&mpy,p); + dbl_mp(e[i].d,&mperr,p); + add(&mpy,&mperr,&mpy1,p); sub(&mpy,&mperr,&mpy2,p); + mp_dbl(&mpy1,&y1,p); mp_dbl(&mpy2,&y2,p); + if (y1==y2) return y1; + } + return y1; } |