]>
Commit | Line | Data |
---|---|---|
f7eac6eb RM |
1 | /* @(#)k_rem_pio2.c 5.1 93/09/24 */ |
2 | /* | |
3 | * ==================================================== | |
4 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | |
5 | * | |
6 | * Developed at SunPro, a Sun Microsystems, Inc. business. | |
7 | * Permission to use, copy, modify, and distribute this | |
6d52618b | 8 | * software is freely granted, provided that this notice |
f7eac6eb RM |
9 | * is preserved. |
10 | * ==================================================== | |
11 | */ | |
12 | ||
13 | #if defined(LIBM_SCCS) && !defined(lint) | |
14 | static char rcsid[] = "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $"; | |
15 | #endif | |
16 | ||
17 | /* | |
18 | * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2) | |
19 | * double x[],y[]; int e0,nx,prec; int ipio2[]; | |
6d52618b UD |
20 | * |
21 | * __kernel_rem_pio2 return the last three digits of N with | |
f7eac6eb RM |
22 | * y = x - N*pi/2 |
23 | * so that |y| < pi/2. | |
24 | * | |
6d52618b | 25 | * The method is to compute the integer (mod 8) and fraction parts of |
f7eac6eb RM |
26 | * (2/pi)*x without doing the full multiplication. In general we |
27 | * skip the part of the product that are known to be a huge integer ( | |
28 | * more accurately, = 0 mod 8 ). Thus the number of operations are | |
29 | * independent of the exponent of the input. | |
30 | * | |
31 | * (2/pi) is represented by an array of 24-bit integers in ipio2[]. | |
32 | * | |
33 | * Input parameters: | |
6d52618b | 34 | * x[] The input value (must be positive) is broken into nx |
f7eac6eb | 35 | * pieces of 24-bit integers in double precision format. |
6d52618b UD |
36 | * x[i] will be the i-th 24 bit of x. The scaled exponent |
37 | * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0 | |
f7eac6eb RM |
38 | * match x's up to 24 bits. |
39 | * | |
40 | * Example of breaking a double positive z into x[0]+x[1]+x[2]: | |
41 | * e0 = ilogb(z)-23 | |
42 | * z = scalbn(z,-e0) | |
43 | * for i = 0,1,2 | |
44 | * x[i] = floor(z) | |
45 | * z = (z-x[i])*2**24 | |
46 | * | |
47 | * | |
48 | * y[] ouput result in an array of double precision numbers. | |
49 | * The dimension of y[] is: | |
50 | * 24-bit precision 1 | |
51 | * 53-bit precision 2 | |
52 | * 64-bit precision 2 | |
53 | * 113-bit precision 3 | |
54 | * The actual value is the sum of them. Thus for 113-bit | |
6d52618b | 55 | * precision, one may have to do something like: |
f7eac6eb RM |
56 | * |
57 | * long double t,w,r_head, r_tail; | |
58 | * t = (long double)y[2] + (long double)y[1]; | |
59 | * w = (long double)y[0]; | |
60 | * r_head = t+w; | |
61 | * r_tail = w - (r_head - t); | |
62 | * | |
63 | * e0 The exponent of x[0] | |
64 | * | |
65 | * nx dimension of x[] | |
66 | * | |
67 | * prec an integer indicating the precision: | |
68 | * 0 24 bits (single) | |
69 | * 1 53 bits (double) | |
70 | * 2 64 bits (extended) | |
71 | * 3 113 bits (quad) | |
72 | * | |
73 | * ipio2[] | |
6d52618b UD |
74 | * integer array, contains the (24*i)-th to (24*i+23)-th |
75 | * bit of 2/pi after binary point. The corresponding | |
f7eac6eb RM |
76 | * floating value is |
77 | * | |
78 | * ipio2[i] * 2^(-24(i+1)). | |
79 | * | |
80 | * External function: | |
81 | * double scalbn(), floor(); | |
82 | * | |
83 | * | |
84 | * Here is the description of some local variables: | |
85 | * | |
86 | * jk jk+1 is the initial number of terms of ipio2[] needed | |
87 | * in the computation. The recommended value is 2,3,4, | |
88 | * 6 for single, double, extended,and quad. | |
89 | * | |
6d52618b UD |
90 | * jz local integer variable indicating the number of |
91 | * terms of ipio2[] used. | |
f7eac6eb RM |
92 | * |
93 | * jx nx - 1 | |
94 | * | |
95 | * jv index for pointing to the suitable ipio2[] for the | |
96 | * computation. In general, we want | |
97 | * ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8 | |
98 | * is an integer. Thus | |
99 | * e0-3-24*jv >= 0 or (e0-3)/24 >= jv | |
100 | * Hence jv = max(0,(e0-3)/24). | |
101 | * | |
102 | * jp jp+1 is the number of terms in PIo2[] needed, jp = jk. | |
103 | * | |
104 | * q[] double array with integral value, representing the | |
105 | * 24-bits chunk of the product of x and 2/pi. | |
106 | * | |
107 | * q0 the corresponding exponent of q[0]. Note that the | |
108 | * exponent for q[i] would be q0-24*i. | |
109 | * | |
110 | * PIo2[] double precision array, obtained by cutting pi/2 | |
6d52618b | 111 | * into 24 bits chunks. |
f7eac6eb | 112 | * |
6d52618b | 113 | * f[] ipio2[] in floating point |
f7eac6eb RM |
114 | * |
115 | * iq[] integer array by breaking up q[] in 24-bits chunk. | |
116 | * | |
117 | * fq[] final product of x*(2/pi) in fq[0],..,fq[jk] | |
118 | * | |
119 | * ih integer. If >0 it indicates q[] is >= 0.5, hence | |
120 | * it also indicates the *sign* of the result. | |
121 | * | |
122 | */ | |
123 | ||
124 | ||
125 | /* | |
126 | * Constants: | |
6d52618b UD |
127 | * The hexadecimal values are the intended ones for the following |
128 | * constants. The decimal values may be used, provided that the | |
129 | * compiler will convert from decimal to binary accurately enough | |
f7eac6eb RM |
130 | * to produce the hexadecimal values shown. |
131 | */ | |
132 | ||
1ed0291c | 133 | #include <math.h> |
aaee3cd8 | 134 | #include <math-narrow-eval.h> |
1ed0291c | 135 | #include <math_private.h> |
9090848d | 136 | #include <libc-diag.h> |
f7eac6eb | 137 | |
f7eac6eb | 138 | static const int init_jk[] = {2,3,4,6}; /* initial value for jk */ |
f7eac6eb | 139 | |
f7eac6eb | 140 | static const double PIo2[] = { |
f7eac6eb RM |
141 | 1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */ |
142 | 7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */ | |
143 | 5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */ | |
144 | 3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */ | |
145 | 1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */ | |
146 | 1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */ | |
147 | 2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */ | |
148 | 2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */ | |
149 | }; | |
150 | ||
6d52618b | 151 | static const double |
c5d5d574 OB |
152 | zero = 0.0, |
153 | one = 1.0, | |
154 | two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ | |
155 | twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */ | |
f7eac6eb | 156 | |
c5d5d574 OB |
157 | int |
158 | __kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec, | |
159 | const int32_t *ipio2) | |
f7eac6eb | 160 | { |
c5d5d574 OB |
161 | int32_t jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih; |
162 | double z, fw, f[20], fq[20], q[20]; | |
f7eac6eb | 163 | |
c5d5d574 OB |
164 | /* initialize jk*/ |
165 | jk = init_jk[prec]; | |
166 | jp = jk; | |
f7eac6eb | 167 | |
c5d5d574 OB |
168 | /* determine jx,jv,q0, note that 3>q0 */ |
169 | jx = nx - 1; | |
170 | jv = (e0 - 3) / 24; if (jv < 0) | |
171 | jv = 0; | |
172 | q0 = e0 - 24 * (jv + 1); | |
f7eac6eb | 173 | |
c5d5d574 OB |
174 | /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */ |
175 | j = jv - jx; m = jx + jk; | |
176 | for (i = 0; i <= m; i++, j++) | |
177 | f[i] = (j < 0) ? zero : (double) ipio2[j]; | |
f7eac6eb | 178 | |
c5d5d574 OB |
179 | /* compute q[0],q[1],...q[jk] */ |
180 | for (i = 0; i <= jk; i++) | |
181 | { | |
182 | for (j = 0, fw = 0.0; j <= jx; j++) | |
183 | fw += x[j] * f[jx + i - j]; | |
184 | q[i] = fw; | |
185 | } | |
f7eac6eb | 186 | |
c5d5d574 | 187 | jz = jk; |
f7eac6eb | 188 | recompute: |
c5d5d574 OB |
189 | /* distill q[] into iq[] reversingly */ |
190 | for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) | |
191 | { | |
192 | fw = (double) ((int32_t) (twon24 * z)); | |
193 | iq[i] = (int32_t) (z - two24 * fw); | |
194 | z = q[j - 1] + fw; | |
195 | } | |
f7eac6eb | 196 | |
c5d5d574 OB |
197 | /* compute n */ |
198 | z = __scalbn (z, q0); /* actual value of z */ | |
e44acb20 | 199 | z -= 8.0 * floor (z * 0.125); /* trim off integer >= 8 */ |
c5d5d574 OB |
200 | n = (int32_t) z; |
201 | z -= (double) n; | |
202 | ih = 0; | |
203 | if (q0 > 0) /* need iq[jz-1] to determine n */ | |
204 | { | |
205 | i = (iq[jz - 1] >> (24 - q0)); n += i; | |
206 | iq[jz - 1] -= i << (24 - q0); | |
207 | ih = iq[jz - 1] >> (23 - q0); | |
208 | } | |
209 | else if (q0 == 0) | |
210 | ih = iq[jz - 1] >> 23; | |
211 | else if (z >= 0.5) | |
212 | ih = 2; | |
f7eac6eb | 213 | |
c5d5d574 OB |
214 | if (ih > 0) /* q > 0.5 */ |
215 | { | |
216 | n += 1; carry = 0; | |
217 | for (i = 0; i < jz; i++) /* compute 1-q */ | |
218 | { | |
219 | j = iq[i]; | |
220 | if (carry == 0) | |
221 | { | |
222 | if (j != 0) | |
223 | { | |
224 | carry = 1; iq[i] = 0x1000000 - j; | |
225 | } | |
f7eac6eb | 226 | } |
c5d5d574 OB |
227 | else |
228 | iq[i] = 0xffffff - j; | |
229 | } | |
230 | if (q0 > 0) /* rare case: chance is 1 in 12 */ | |
231 | { | |
232 | switch (q0) | |
233 | { | |
234 | case 1: | |
235 | iq[jz - 1] &= 0x7fffff; break; | |
236 | case 2: | |
237 | iq[jz - 1] &= 0x3fffff; break; | |
f7eac6eb RM |
238 | } |
239 | } | |
c5d5d574 OB |
240 | if (ih == 2) |
241 | { | |
242 | z = one - z; | |
243 | if (carry != 0) | |
244 | z -= __scalbn (one, q0); | |
245 | } | |
246 | } | |
f7eac6eb | 247 | |
c5d5d574 OB |
248 | /* check if recomputation is needed */ |
249 | if (z == zero) | |
250 | { | |
251 | j = 0; | |
252 | for (i = jz - 1; i >= jk; i--) | |
253 | j |= iq[i]; | |
254 | if (j == 0) /* need recomputation */ | |
255 | { | |
b65f0b7b SL |
256 | /* On s390x gcc 6.1 -O3 produces the warning "array subscript is below |
257 | array bounds [-Werror=array-bounds]". Only __ieee754_rem_pio2l | |
258 | calls __kernel_rem_pio2 for normal numbers and |x| > pi/4 in case | |
259 | of ldbl-96 and |x| > 3pi/4 in case of ldbl-128[ibm]. | |
260 | Thus x can't be zero and ipio2 is not zero, too. Thus not all iq[] | |
261 | values can't be zero. */ | |
262 | DIAG_PUSH_NEEDS_COMMENT; | |
263 | DIAG_IGNORE_NEEDS_COMMENT (6.1, "-Warray-bounds"); | |
c5d5d574 OB |
264 | for (k = 1; iq[jk - k] == 0; k++) |
265 | ; /* k = no. of terms needed */ | |
b65f0b7b | 266 | DIAG_POP_NEEDS_COMMENT; |
f7eac6eb | 267 | |
c5d5d574 OB |
268 | for (i = jz + 1; i <= jz + k; i++) /* add q[jz+1] to q[jz+k] */ |
269 | { | |
270 | f[jx + i] = (double) ipio2[jv + i]; | |
271 | for (j = 0, fw = 0.0; j <= jx; j++) | |
272 | fw += x[j] * f[jx + i - j]; | |
273 | q[i] = fw; | |
f7eac6eb | 274 | } |
c5d5d574 OB |
275 | jz += k; |
276 | goto recompute; | |
f7eac6eb | 277 | } |
c5d5d574 | 278 | } |
f7eac6eb | 279 | |
c5d5d574 OB |
280 | /* chop off zero terms */ |
281 | if (z == 0.0) | |
282 | { | |
283 | jz -= 1; q0 -= 24; | |
284 | while (iq[jz] == 0) | |
285 | { | |
286 | jz--; q0 -= 24; | |
f7eac6eb | 287 | } |
c5d5d574 OB |
288 | } |
289 | else /* break z into 24-bit if necessary */ | |
290 | { | |
291 | z = __scalbn (z, -q0); | |
292 | if (z >= two24) | |
293 | { | |
294 | fw = (double) ((int32_t) (twon24 * z)); | |
295 | iq[jz] = (int32_t) (z - two24 * fw); | |
296 | jz += 1; q0 += 24; | |
297 | iq[jz] = (int32_t) fw; | |
f7eac6eb | 298 | } |
c5d5d574 OB |
299 | else |
300 | iq[jz] = (int32_t) z; | |
301 | } | |
f7eac6eb | 302 | |
c5d5d574 OB |
303 | /* convert integer "bit" chunk to floating-point value */ |
304 | fw = __scalbn (one, q0); | |
305 | for (i = jz; i >= 0; i--) | |
306 | { | |
307 | q[i] = fw * (double) iq[i]; fw *= twon24; | |
308 | } | |
f7eac6eb | 309 | |
c5d5d574 OB |
310 | /* compute PIo2[0,...,jp]*q[jz,...,0] */ |
311 | for (i = jz; i >= 0; i--) | |
312 | { | |
313 | for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++) | |
314 | fw += PIo2[k] * q[i + k]; | |
315 | fq[jz - i] = fw; | |
316 | } | |
317 | ||
318 | /* compress fq[] into y[] */ | |
319 | switch (prec) | |
320 | { | |
321 | case 0: | |
322 | fw = 0.0; | |
323 | for (i = jz; i >= 0; i--) | |
324 | fw += fq[i]; | |
325 | y[0] = (ih == 0) ? fw : -fw; | |
326 | break; | |
327 | case 1: | |
328 | case 2:; | |
c5d5d574 OB |
329 | double fv = 0.0; |
330 | for (i = jz; i >= 0; i--) | |
54142c44 | 331 | fv = math_narrow_eval (fv + fq[i]); |
c5d5d574 | 332 | y[0] = (ih == 0) ? fv : -fv; |
3d6302a5 JM |
333 | /* GCC mainline (to be GCC 9), as of 2018-05-22 on i686, warns |
334 | that fq[0] may be used uninitialized. This is not possible | |
335 | because jz is always nonnegative when the above loop | |
336 | initializing fq is executed, because the result is never zero | |
337 | to full precision (this function is not called for zero | |
338 | arguments). */ | |
339 | DIAG_PUSH_NEEDS_COMMENT; | |
340 | DIAG_IGNORE_NEEDS_COMMENT (9, "-Wmaybe-uninitialized"); | |
54142c44 | 341 | fv = math_narrow_eval (fq[0] - fv); |
3d6302a5 | 342 | DIAG_POP_NEEDS_COMMENT; |
c5d5d574 | 343 | for (i = 1; i <= jz; i++) |
54142c44 | 344 | fv = math_narrow_eval (fv + fq[i]); |
c5d5d574 OB |
345 | y[1] = (ih == 0) ? fv : -fv; |
346 | break; | |
347 | case 3: /* painful */ | |
348 | for (i = jz; i > 0; i--) | |
349 | { | |
54142c44 | 350 | double fv = math_narrow_eval (fq[i - 1] + fq[i]); |
c5d5d574 OB |
351 | fq[i] += fq[i - 1] - fv; |
352 | fq[i - 1] = fv; | |
353 | } | |
354 | for (i = jz; i > 1; i--) | |
355 | { | |
54142c44 | 356 | double fv = math_narrow_eval (fq[i - 1] + fq[i]); |
c5d5d574 OB |
357 | fq[i] += fq[i - 1] - fv; |
358 | fq[i - 1] = fv; | |
359 | } | |
360 | for (fw = 0.0, i = jz; i >= 2; i--) | |
361 | fw += fq[i]; | |
362 | if (ih == 0) | |
363 | { | |
364 | y[0] = fq[0]; y[1] = fq[1]; y[2] = fw; | |
365 | } | |
366 | else | |
367 | { | |
368 | y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw; | |
f7eac6eb | 369 | } |
c5d5d574 OB |
370 | } |
371 | return n & 7; | |
f7eac6eb | 372 | } |