]>
Commit | Line | Data |
---|---|---|
f7eac6eb | 1 | /* |
e4d82761 | 2 | * IBM Accurate Mathematical Library |
aeb25823 | 3 | * written by International Business Machines Corp. |
8ec250a4 | 4 | * Copyright (C) 2001, 2011 Free Software Foundation |
f7eac6eb | 5 | * |
e4d82761 UD |
6 | * This program is free software; you can redistribute it and/or modify |
7 | * it under the terms of the GNU Lesser General Public License as published by | |
cc7375ce | 8 | * the Free Software Foundation; either version 2.1 of the License, or |
e4d82761 | 9 | * (at your option) any later version. |
6d52618b | 10 | * |
e4d82761 UD |
11 | * This program is distributed in the hope that it will be useful, |
12 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | |
13 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
c6c6dd48 | 14 | * GNU Lesser General Public License for more details. |
f7eac6eb | 15 | * |
e4d82761 UD |
16 | * You should have received a copy of the GNU Lesser General Public License |
17 | * along with this program; if not, write to the Free Software | |
18 | * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. | |
f7eac6eb | 19 | */ |
e4d82761 UD |
20 | /*********************************************************************/ |
21 | /* */ | |
aeb25823 | 22 | /* MODULE_NAME:ulog.c */ |
e4d82761 UD |
23 | /* */ |
24 | /* FUNCTION:ulog */ | |
25 | /* */ | |
26 | /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h ulog.h */ | |
27 | /* mpexp.c mplog.c mpa.c */ | |
28 | /* ulog.tbl */ | |
29 | /* */ | |
30 | /* An ultimate log routine. Given an IEEE double machine number x */ | |
31 | /* it computes the correctly rounded (to nearest) value of log(x). */ | |
32 | /* Assumption: Machine arithmetic operations are performed in */ | |
33 | /* round to nearest mode of IEEE 754 standard. */ | |
34 | /* */ | |
35 | /*********************************************************************/ | |
36 | ||
37 | ||
38 | #include "endian.h" | |
c8b3296b | 39 | #include <dla.h> |
e4d82761 UD |
40 | #include "mpa.h" |
41 | #include "MathLib.h" | |
e859d1d9 AJ |
42 | #include "math_private.h" |
43 | ||
e4d82761 UD |
44 | void __mplog(mp_no *, mp_no *, int); |
45 | ||
46 | /*********************************************************************/ | |
47 | /* An ultimate log routine. Given an IEEE double machine number x */ | |
48 | /* it computes the correctly rounded (to nearest) value of log(x). */ | |
49 | /*********************************************************************/ | |
50 | double __ieee754_log(double x) { | |
51 | #define M 4 | |
52 | static const int pr[M]={8,10,18,32}; | |
50944bca UD |
53 | int i,j,n,ux,dx,p; |
54 | #if 0 | |
55 | int k; | |
56 | #endif | |
e4d82761 | 57 | double dbl_n,u,p0,q,r0,w,nln2a,luai,lubi,lvaj,lvbj, |
0ac5ae23 | 58 | sij,ssij,ttij,A,B,B0,y,y1,y2,polI,polII,sa,sb, |
a1a87169 | 59 | t1,t2,t7,t8,t,ra,rb,ww, |
0ac5ae23 | 60 | a0,aa0,s1,s2,ss2,s3,ss3,a1,aa1,a,aa,b,bb,c; |
a1a87169 UD |
61 | #ifndef DLA_FMA |
62 | double t3,t4,t5,t6; | |
63 | #endif | |
e4d82761 UD |
64 | number num; |
65 | mp_no mpx,mpy,mpy1,mpy2,mperr; | |
66 | ||
67 | #include "ulog.tbl" | |
68 | #include "ulog.h" | |
69 | ||
70 | /* Treating special values of x ( x<=0, x=INF, x=NaN etc.). */ | |
71 | ||
72 | num.d = x; ux = num.i[HIGH_HALF]; dx = num.i[LOW_HALF]; | |
73 | n=0; | |
8ec250a4 | 74 | if (__builtin_expect(ux < 0x00100000, 0)) { |
0ac5ae23 UD |
75 | if (__builtin_expect(((ux & 0x7fffffff) | dx) == 0, 0)) |
76 | return MHALF/ZERO; /* return -INF */ | |
77 | if (__builtin_expect(ux < 0, 0)) | |
78 | return (x-x)/ZERO; /* return NaN */ | |
e4d82761 UD |
79 | n -= 54; x *= two54.d; /* scale x */ |
80 | num.d = x; | |
81 | } | |
0ac5ae23 UD |
82 | if (__builtin_expect(ux >= 0x7ff00000, 0)) |
83 | return x+x; /* INF or NaN */ | |
e4d82761 UD |
84 | |
85 | /* Regular values of x */ | |
86 | ||
87 | w = x-ONE; | |
8ec250a4 | 88 | if (__builtin_expect(ABS(w) > U03, 1)) { goto case_03; } |
e4d82761 UD |
89 | |
90 | ||
91 | /*--- Stage I, the case abs(x-1) < 0.03 */ | |
92 | ||
93 | t8 = MHALF*w; | |
94 | EMULV(t8,w,a,aa,t1,t2,t3,t4,t5) | |
95 | EADD(w,a,b,bb) | |
96 | ||
97 | /* Evaluate polynomial II */ | |
98 | polII = (b0.d+w*(b1.d+w*(b2.d+w*(b3.d+w*(b4.d+ | |
0ac5ae23 | 99 | w*(b5.d+w*(b6.d+w*(b7.d+w*b8.d))))))))*w*w*w; |
e4d82761 UD |
100 | c = (aa+bb)+polII; |
101 | ||
102 | /* End stage I, case abs(x-1) < 0.03 */ | |
103 | if ((y=b+(c+b*E2)) == b+(c-b*E2)) return y; | |
104 | ||
105 | /*--- Stage II, the case abs(x-1) < 0.03 */ | |
106 | ||
107 | a = d11.d+w*(d12.d+w*(d13.d+w*(d14.d+w*(d15.d+w*(d16.d+ | |
0ac5ae23 | 108 | w*(d17.d+w*(d18.d+w*(d19.d+w*d20.d)))))))); |
e4d82761 UD |
109 | EMULV(w,a,s2,ss2,t1,t2,t3,t4,t5) |
110 | ADD2(d10.d,dd10.d,s2,ss2,s3,ss3,t1,t2) | |
111 | MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) | |
112 | ADD2(d9.d,dd9.d,s2,ss2,s3,ss3,t1,t2) | |
113 | MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) | |
114 | ADD2(d8.d,dd8.d,s2,ss2,s3,ss3,t1,t2) | |
115 | MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) | |
116 | ADD2(d7.d,dd7.d,s2,ss2,s3,ss3,t1,t2) | |
117 | MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) | |
118 | ADD2(d6.d,dd6.d,s2,ss2,s3,ss3,t1,t2) | |
119 | MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) | |
120 | ADD2(d5.d,dd5.d,s2,ss2,s3,ss3,t1,t2) | |
121 | MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) | |
122 | ADD2(d4.d,dd4.d,s2,ss2,s3,ss3,t1,t2) | |
123 | MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) | |
124 | ADD2(d3.d,dd3.d,s2,ss2,s3,ss3,t1,t2) | |
125 | MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) | |
126 | ADD2(d2.d,dd2.d,s2,ss2,s3,ss3,t1,t2) | |
127 | MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) | |
128 | MUL2(w,ZERO,s2,ss2,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8) | |
129 | ADD2(w,ZERO, s3,ss3, b, bb,t1,t2) | |
130 | ||
131 | /* End stage II, case abs(x-1) < 0.03 */ | |
132 | if ((y=b+(bb+b*E4)) == b+(bb-b*E4)) return y; | |
133 | goto stage_n; | |
134 | ||
135 | /*--- Stage I, the case abs(x-1) > 0.03 */ | |
136 | case_03: | |
137 | ||
138 | /* Find n,u such that x = u*2**n, 1/sqrt(2) < u < sqrt(2) */ | |
139 | n += (num.i[HIGH_HALF] >> 20) - 1023; | |
140 | num.i[HIGH_HALF] = (num.i[HIGH_HALF] & 0x000fffff) | 0x3ff00000; | |
141 | if (num.d > SQRT_2) { num.d *= HALF; n++; } | |
142 | u = num.d; dbl_n = (double) n; | |
143 | ||
144 | /* Find i such that ui=1+(i-75)/2**8 is closest to u (i= 0,1,2,...,181) */ | |
145 | num.d += h1.d; | |
146 | i = (num.i[HIGH_HALF] & 0x000fffff) >> 12; | |
147 | ||
148 | /* Find j such that vj=1+(j-180)/2**16 is closest to v=u/ui (j= 0,...,361) */ | |
149 | num.d = u*Iu[i].d + h2.d; | |
150 | j = (num.i[HIGH_HALF] & 0x000fffff) >> 4; | |
151 | ||
152 | /* Compute w=(u-ui*vj)/(ui*vj) */ | |
153 | p0=(ONE+(i-75)*DEL_U)*(ONE+(j-180)*DEL_V); | |
154 | q=u-p0; r0=Iu[i].d*Iv[j].d; w=q*r0; | |
155 | ||
156 | /* Evaluate polynomial I */ | |
157 | polI = w+(a2.d+a3.d*w)*w*w; | |
158 | ||
159 | /* Add up everything */ | |
160 | nln2a = dbl_n*LN2A; | |
161 | luai = Lu[i][0].d; lubi = Lu[i][1].d; | |
162 | lvaj = Lv[j][0].d; lvbj = Lv[j][1].d; | |
163 | EADD(luai,lvaj,sij,ssij) | |
164 | EADD(nln2a,sij,A ,ttij) | |
165 | B0 = (((lubi+lvbj)+ssij)+ttij)+dbl_n*LN2B; | |
166 | B = polI+B0; | |
167 | ||
168 | /* End stage I, case abs(x-1) >= 0.03 */ | |
169 | if ((y=A+(B+E1)) == A+(B-E1)) return y; | |
170 | ||
171 | ||
172 | /*--- Stage II, the case abs(x-1) > 0.03 */ | |
173 | ||
174 | /* Improve the accuracy of r0 */ | |
175 | EMULV(p0,r0,sa,sb,t1,t2,t3,t4,t5) | |
176 | t=r0*((ONE-sa)-sb); | |
177 | EADD(r0,t,ra,rb) | |
178 | ||
179 | /* Compute w */ | |
180 | MUL2(q,ZERO,ra,rb,w,ww,t1,t2,t3,t4,t5,t6,t7,t8) | |
181 | ||
182 | EADD(A,B0,a0,aa0) | |
183 | ||
184 | /* Evaluate polynomial III */ | |
185 | s1 = (c3.d+(c4.d+c5.d*w)*w)*w; | |
186 | EADD(c2.d,s1,s2,ss2) | |
187 | MUL2(s2,ss2,w,ww,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8) | |
188 | MUL2(s3,ss3,w,ww,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) | |
189 | ADD2(s2,ss2,w,ww,s3,ss3,t1,t2) | |
190 | ADD2(s3,ss3,a0,aa0,a1,aa1,t1,t2) | |
191 | ||
192 | /* End stage II, case abs(x-1) >= 0.03 */ | |
193 | if ((y=a1+(aa1+E3)) == a1+(aa1-E3)) return y; | |
194 | ||
195 | ||
196 | /* Final stages. Use multi-precision arithmetic. */ | |
197 | stage_n: | |
f7eac6eb | 198 | |
e4d82761 UD |
199 | for (i=0; i<M; i++) { |
200 | p = pr[i]; | |
ca58f1db | 201 | __dbl_mp(x,&mpx,p); __dbl_mp(y,&mpy,p); |
e4d82761 | 202 | __mplog(&mpx,&mpy,p); |
ca58f1db UD |
203 | __dbl_mp(e[i].d,&mperr,p); |
204 | __add(&mpy,&mperr,&mpy1,p); __sub(&mpy,&mperr,&mpy2,p); | |
50944bca | 205 | __mp_dbl(&mpy1,&y1,p); __mp_dbl(&mpy2,&y2,p); |
e4d82761 UD |
206 | if (y1==y2) return y1; |
207 | } | |
208 | return y1; | |
f7eac6eb | 209 | } |
0ac5ae23 | 210 | strong_alias (__ieee754_log, __log_finite) |