1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2018 Free Software Foundation, Inc.
5 *
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.
10 *
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.
15 *
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/>.
18 */
19/*********************************************************************/
20/* */
21/* MODULE_NAME:ulog.c */
22/* */
23/* FUNCTION:ulog */
24/* */
25/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h ulog.h */
26/* ulog.tbl */
27/* */
28/* An ultimate log routine. Given an IEEE double machine number x */
29/* it computes the rounded (to nearest) value of log(x). */
30/* Assumption: Machine arithmetic operations are performed in */
31/* round to nearest mode of IEEE 754 standard. */
32/* */
33/*********************************************************************/
34
35
36#include "endian.h"
37#include <dla.h>
38#include "mpa.h"
39#include "MathLib.h"
40#include <math.h>
41#include <math_private.h>
42
43#ifndef SECTION
44# define SECTION
45#endif
46
47/*********************************************************************/
48/* An ultimate log routine. Given an IEEE double machine number x */
49/* it computes the rounded (to nearest) value of log(x). */
50/*********************************************************************/
51double
52SECTION
53__ieee754_log (double x)
54{
55 int i, j, n, ux, dx;
56 double dbl_n, u, p0, q, r0, w, nln2a, luai, lubi, lvaj, lvbj,
57 sij, ssij, ttij, A, B, B0, polI, polII, t8, a, aa, b, bb, c;
58#ifndef DLA_FMS
59 double t1, t2, t3, t4, t5;
60#endif
61 number num;
62
63#include "ulog.tbl"
64#include "ulog.h"
65
66 /* Treating special values of x ( x<=0, x=INF, x=NaN etc.). */
67
68 num.d = x;
69 ux = num.i[HIGH_HALF];
70 dx = num.i[LOW_HALF];
71 n = 0;
72 if (__glibc_unlikely (ux < 0x00100000))
73 {
74 if (__glibc_unlikely (((ux & 0x7fffffff) | dx) == 0))
75 return MHALF / 0.0; /* return -INF */
76 if (__glibc_unlikely (ux < 0))
77 return (x - x) / 0.0; /* return NaN */
78 n -= 54;
79 x *= two54.d; /* scale x */
80 num.d = x;
81 }
82 if (__glibc_unlikely (ux >= 0x7ff00000))
83 return x + x; /* INF or NaN */
84
85 /* Regular values of x */
86
87 w = x - 1;
88 if (__glibc_likely (fabs (w) > U03))
89 goto case_03;
90
91 /* log (1) is +0 in all rounding modes. */
92 if (w == 0.0)
93 return 0.0;
94
95 /*--- The case abs(x-1) < 0.03 */
96
97 t8 = MHALF * w;
98 EMULV (t8, w, a, aa, t1, t2, t3, t4, t5);
99 EADD (w, a, b, bb);
100 /* Evaluate polynomial II */
101 polII = b7.d + w * b8.d;
102 polII = b6.d + w * polII;
103 polII = b5.d + w * polII;
104 polII = b4.d + w * polII;
105 polII = b3.d + w * polII;
106 polII = b2.d + w * polII;
107 polII = b1.d + w * polII;
108 polII = b0.d + w * polII;
109 polII *= w * w * w;
110 c = (aa + bb) + polII;
111
112 /* Here b contains the high part of the result, and c the low part.
113 Maximum error is b * 2.334e-19, so accuracy is >61 bits.
114 Therefore max ULP error of b + c is ~0.502. */
115 return b + c;
116
117 /*--- The case abs(x-1) > 0.03 */
118case_03:
119
120 /* Find n,u such that x = u*2**n, 1/sqrt(2) < u < sqrt(2) */
121 n += (num.i[HIGH_HALF] >> 20) - 1023;
122 num.i[HIGH_HALF] = (num.i[HIGH_HALF] & 0x000fffff) | 0x3ff00000;
123 if (num.d > SQRT_2)
124 {
125 num.d *= HALF;
126 n++;
127 }
128 u = num.d;
129 dbl_n = (double) n;
130
131 /* Find i such that ui=1+(i-75)/2**8 is closest to u (i= 0,1,2,...,181) */
132 num.d += h1.d;
133 i = (num.i[HIGH_HALF] & 0x000fffff) >> 12;
134
135 /* Find j such that vj=1+(j-180)/2**16 is closest to v=u/ui (j= 0,...,361) */
136 num.d = u * Iu[i].d + h2.d;
137 j = (num.i[HIGH_HALF] & 0x000fffff) >> 4;
138
139 /* Compute w=(u-ui*vj)/(ui*vj) */
140 p0 = (1 + (i - 75) * DEL_U) * (1 + (j - 180) * DEL_V);
141 q = u - p0;
142 r0 = Iu[i].d * Iv[j].d;
143 w = q * r0;
144
145 /* Evaluate polynomial I */
146 polI = w + (a2.d + a3.d * w) * w * w;
147
148 /* Add up everything */
149 nln2a = dbl_n * LN2A;
150 luai = Lu[i][0].d;
151 lubi = Lu[i][1].d;
152 lvaj = Lv[j][0].d;
153 lvbj = Lv[j][1].d;
154 EADD (luai, lvaj, sij, ssij);
155 EADD (nln2a, sij, A, ttij);
156 B0 = (((lubi + lvbj) + ssij) + ttij) + dbl_n * LN2B;
157 B = polI + B0;
158
159 /* Here A contains the high part of the result, and B the low part.
160 Maximum abs error is 6.095e-21 and min log (x) is 0.0295 since x > 1.03.
161 Therefore max ULP error of A + B is ~0.502. */
162 return A + B;
163}
164
165#ifndef __ieee754_log
166strong_alias (__ieee754_log, __log_finite)
167#endif
168