2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2014 Free Software Foundation, Inc.
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
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
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
14 * GNU Lesser General Public License for more details.
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, see <http://www.gnu.org/licenses/>.
19 /*********************************************************************/
21 /* MODULE_NAME:ulog.c */
25 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h ulog.h */
26 /* mpexp.c mplog.c mpa.c */
29 /* An ultimate log routine. Given an IEEE double machine number x */
30 /* it computes the correctly rounded (to nearest) value of log(x). */
31 /* Assumption: Machine arithmetic operations are performed in */
32 /* round to nearest mode of IEEE 754 standard. */
34 /*********************************************************************/
41 #include <math_private.h>
42 #include <stap-probe.h>
48 void __mplog (mp_no
*, mp_no
*, int);
50 /*********************************************************************/
51 /* An ultimate log routine. Given an IEEE double machine number x */
52 /* it computes the correctly rounded (to nearest) value of log(x). */
53 /*********************************************************************/
56 __ieee754_log (double x
)
59 static const int pr
[M
] = { 8, 10, 18, 32 };
60 int i
, j
, n
, ux
, dx
, p
;
61 double dbl_n
, u
, p0
, q
, r0
, w
, nln2a
, luai
, lubi
, lvaj
, lvbj
,
62 sij
, ssij
, ttij
, A
, B
, B0
, y
, y1
, y2
, polI
, polII
, sa
, sb
,
63 t1
, t2
, t7
, t8
, t
, ra
, rb
, ww
,
64 a0
, aa0
, s1
, s2
, ss2
, s3
, ss3
, a1
, aa1
, a
, aa
, b
, bb
, c
;
66 double t3
, t4
, t5
, t6
;
69 mp_no mpx
, mpy
, mpy1
, mpy2
, mperr
;
74 /* Treating special values of x ( x<=0, x=INF, x=NaN etc.). */
77 ux
= num
.i
[HIGH_HALF
];
80 if (__builtin_expect (ux
< 0x00100000, 0))
82 if (__builtin_expect (((ux
& 0x7fffffff) | dx
) == 0, 0))
83 return MHALF
/ 0.0; /* return -INF */
84 if (__builtin_expect (ux
< 0, 0))
85 return (x
- x
) / 0.0; /* return NaN */
87 x
*= two54
.d
; /* scale x */
90 if (__builtin_expect (ux
>= 0x7ff00000, 0))
91 return x
+ x
; /* INF or NaN */
93 /* Regular values of x */
96 if (__builtin_expect (ABS (w
) > U03
, 1))
99 /*--- Stage I, the case abs(x-1) < 0.03 */
102 EMULV (t8
, w
, a
, aa
, t1
, t2
, t3
, t4
, t5
);
104 /* Evaluate polynomial II */
105 polII
= b7
.d
+ w
* b8
.d
;
106 polII
= b6
.d
+ w
* polII
;
107 polII
= b5
.d
+ w
* polII
;
108 polII
= b4
.d
+ w
* polII
;
109 polII
= b3
.d
+ w
* polII
;
110 polII
= b2
.d
+ w
* polII
;
111 polII
= b1
.d
+ w
* polII
;
112 polII
= b0
.d
+ w
* polII
;
114 c
= (aa
+ bb
) + polII
;
116 /* End stage I, case abs(x-1) < 0.03 */
117 if ((y
= b
+ (c
+ b
* E2
)) == b
+ (c
- b
* E2
))
120 /*--- Stage II, the case abs(x-1) < 0.03 */
122 a
= d19
.d
+ w
* d20
.d
;
132 EMULV (w
, a
, s2
, ss2
, t1
, t2
, t3
, t4
, t5
);
133 ADD2 (d10
.d
, dd10
.d
, s2
, ss2
, s3
, ss3
, t1
, t2
);
134 MUL2 (w
, 0, s3
, ss3
, s2
, ss2
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
135 ADD2 (d9
.d
, dd9
.d
, s2
, ss2
, s3
, ss3
, t1
, t2
);
136 MUL2 (w
, 0, s3
, ss3
, s2
, ss2
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
137 ADD2 (d8
.d
, dd8
.d
, s2
, ss2
, s3
, ss3
, t1
, t2
);
138 MUL2 (w
, 0, s3
, ss3
, s2
, ss2
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
139 ADD2 (d7
.d
, dd7
.d
, s2
, ss2
, s3
, ss3
, t1
, t2
);
140 MUL2 (w
, 0, s3
, ss3
, s2
, ss2
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
141 ADD2 (d6
.d
, dd6
.d
, s2
, ss2
, s3
, ss3
, t1
, t2
);
142 MUL2 (w
, 0, s3
, ss3
, s2
, ss2
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
143 ADD2 (d5
.d
, dd5
.d
, s2
, ss2
, s3
, ss3
, t1
, t2
);
144 MUL2 (w
, 0, s3
, ss3
, s2
, ss2
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
145 ADD2 (d4
.d
, dd4
.d
, s2
, ss2
, s3
, ss3
, t1
, t2
);
146 MUL2 (w
, 0, s3
, ss3
, s2
, ss2
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
147 ADD2 (d3
.d
, dd3
.d
, s2
, ss2
, s3
, ss3
, t1
, t2
);
148 MUL2 (w
, 0, s3
, ss3
, s2
, ss2
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
149 ADD2 (d2
.d
, dd2
.d
, s2
, ss2
, s3
, ss3
, t1
, t2
);
150 MUL2 (w
, 0, s3
, ss3
, s2
, ss2
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
151 MUL2 (w
, 0, s2
, ss2
, s3
, ss3
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
152 ADD2 (w
, 0, s3
, ss3
, b
, bb
, t1
, t2
);
154 /* End stage II, case abs(x-1) < 0.03 */
155 if ((y
= b
+ (bb
+ b
* E4
)) == b
+ (bb
- b
* E4
))
159 /*--- Stage I, the case abs(x-1) > 0.03 */
162 /* Find n,u such that x = u*2**n, 1/sqrt(2) < u < sqrt(2) */
163 n
+= (num
.i
[HIGH_HALF
] >> 20) - 1023;
164 num
.i
[HIGH_HALF
] = (num
.i
[HIGH_HALF
] & 0x000fffff) | 0x3ff00000;
173 /* Find i such that ui=1+(i-75)/2**8 is closest to u (i= 0,1,2,...,181) */
175 i
= (num
.i
[HIGH_HALF
] & 0x000fffff) >> 12;
177 /* Find j such that vj=1+(j-180)/2**16 is closest to v=u/ui (j= 0,...,361) */
178 num
.d
= u
* Iu
[i
].d
+ h2
.d
;
179 j
= (num
.i
[HIGH_HALF
] & 0x000fffff) >> 4;
181 /* Compute w=(u-ui*vj)/(ui*vj) */
182 p0
= (1 + (i
- 75) * DEL_U
) * (1 + (j
- 180) * DEL_V
);
184 r0
= Iu
[i
].d
* Iv
[j
].d
;
187 /* Evaluate polynomial I */
188 polI
= w
+ (a2
.d
+ a3
.d
* w
) * w
* w
;
190 /* Add up everything */
191 nln2a
= dbl_n
* LN2A
;
196 EADD (luai
, lvaj
, sij
, ssij
);
197 EADD (nln2a
, sij
, A
, ttij
);
198 B0
= (((lubi
+ lvbj
) + ssij
) + ttij
) + dbl_n
* LN2B
;
201 /* End stage I, case abs(x-1) >= 0.03 */
202 if ((y
= A
+ (B
+ E1
)) == A
+ (B
- E1
))
206 /*--- Stage II, the case abs(x-1) > 0.03 */
208 /* Improve the accuracy of r0 */
209 EMULV (p0
, r0
, sa
, sb
, t1
, t2
, t3
, t4
, t5
);
210 t
= r0
* ((1 - sa
) - sb
);
211 EADD (r0
, t
, ra
, rb
);
214 MUL2 (q
, 0, ra
, rb
, w
, ww
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
216 EADD (A
, B0
, a0
, aa0
);
218 /* Evaluate polynomial III */
219 s1
= (c3
.d
+ (c4
.d
+ c5
.d
* w
) * w
) * w
;
220 EADD (c2
.d
, s1
, s2
, ss2
);
221 MUL2 (s2
, ss2
, w
, ww
, s3
, ss3
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
222 MUL2 (s3
, ss3
, w
, ww
, s2
, ss2
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
223 ADD2 (s2
, ss2
, w
, ww
, s3
, ss3
, t1
, t2
);
224 ADD2 (s3
, ss3
, a0
, aa0
, a1
, aa1
, t1
, t2
);
226 /* End stage II, case abs(x-1) >= 0.03 */
227 if ((y
= a1
+ (aa1
+ E3
)) == a1
+ (aa1
- E3
))
231 /* Final stages. Use multi-precision arithmetic. */
234 for (i
= 0; i
< M
; i
++)
237 __dbl_mp (x
, &mpx
, p
);
238 __dbl_mp (y
, &mpy
, p
);
239 __mplog (&mpx
, &mpy
, p
);
240 __dbl_mp (e
[i
].d
, &mperr
, p
);
241 __add (&mpy
, &mperr
, &mpy1
, p
);
242 __sub (&mpy
, &mperr
, &mpy2
, p
);
243 __mp_dbl (&mpy1
, &y1
, p
);
244 __mp_dbl (&mpy2
, &y2
, p
);
247 LIBC_PROBE (slowlog
, 3, &p
, &x
, &y1
);
251 LIBC_PROBE (slowlog_inexact
, 3, &p
, &x
, &y1
);
255 #ifndef __ieee754_log
256 strong_alias (__ieee754_log
, __log_finite
)