| blob | 958f28c2557caf23ce9476771f1787cf071d19de |
1 /* origin: FreeBSD /usr/src/lib/msun/src/k_rem_pio2.c */
2 /*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunSoft, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12 /*
13 * __rem_pio2_large(x,y,e0,nx,prec)
14 * double x[],y[]; int e0,nx,prec;
15 *
16 * __rem_pio2_large return the last three digits of N with
17 * y = x - N*pi/2
18 * so that |y| < pi/2.
19 *
20 * The method is to compute the integer (mod 8) and fraction parts of
21 * (2/pi)*x without doing the full multiplication. In general we
22 * skip the part of the product that are known to be a huge integer (
23 * more accurately, = 0 mod 8 ). Thus the number of operations are
24 * independent of the exponent of the input.
25 *
26 * (2/pi) is represented by an array of 24-bit integers in ipio2[].
27 *
28 * Input parameters:
29 * x[] The input value (must be positive) is broken into nx
30 * pieces of 24-bit integers in double precision format.
31 * x[i] will be the i-th 24 bit of x. The scaled exponent
32 * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
33 * match x's up to 24 bits.
34 *
35 * Example of breaking a double positive z into x[0]+x[1]+x[2]:
36 * e0 = ilogb(z)-23
37 * z = scalbn(z,-e0)
38 * for i = 0,1,2
39 * x[i] = floor(z)
40 * z = (z-x[i])*2**24
41 *
42 *
43 * y[] ouput result in an array of double precision numbers.
44 * The dimension of y[] is:
45 * 24-bit precision 1
46 * 53-bit precision 2
47 * 64-bit precision 2
48 * 113-bit precision 3
49 * The actual value is the sum of them. Thus for 113-bit
50 * precison, one may have to do something like:
51 *
52 * long double t,w,r_head, r_tail;
53 * t = (long double)y[2] + (long double)y[1];
54 * w = (long double)y[0];
55 * r_head = t+w;
56 * r_tail = w - (r_head - t);
57 *
58 * e0 The exponent of x[0]. Must be <= 16360 or you need to
59 * expand the ipio2 table.
60 *
61 * nx dimension of x[]
62 *
63 * prec an integer indicating the precision:
64 * 0 24 bits (single)
65 * 1 53 bits (double)
66 * 2 64 bits (extended)
67 * 3 113 bits (quad)
68 *
69 * External function:
70 * double scalbn(), floor();
71 *
72 *
73 * Here is the description of some local variables:
74 *
75 * jk jk+1 is the initial number of terms of ipio2[] needed
76 * in the computation. The minimum and recommended value
77 * for jk is 3,4,4,6 for single, double, extended, and quad.
78 * jk+1 must be 2 larger than you might expect so that our
79 * recomputation test works. (Up to 24 bits in the integer
80 * part (the 24 bits of it that we compute) and 23 bits in
81 * the fraction part may be lost to cancelation before we
82 * recompute.)
83 *
84 * jz local integer variable indicating the number of
85 * terms of ipio2[] used.
86 *
87 * jx nx - 1
88 *
89 * jv index for pointing to the suitable ipio2[] for the
90 * computation. In general, we want
91 * ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
92 * is an integer. Thus
93 * e0-3-24*jv >= 0 or (e0-3)/24 >= jv
94 * Hence jv = max(0,(e0-3)/24).
95 *
96 * jp jp+1 is the number of terms in PIo2[] needed, jp = jk.
97 *
98 * q[] double array with integral value, representing the
99 * 24-bits chunk of the product of x and 2/pi.
100 *
101 * q0 the corresponding exponent of q[0]. Note that the
102 * exponent for q[i] would be q0-24*i.
103 *
104 * PIo2[] double precision array, obtained by cutting pi/2
105 * into 24 bits chunks.
106 *
107 * f[] ipio2[] in floating point
108 *
109 * iq[] integer array by breaking up q[] in 24-bits chunk.
110 *
111 * fq[] final product of x*(2/pi) in fq[0],..,fq[jk]
112 *
113 * ih integer. If >0 it indicates q[] is >= 0.5, hence
114 * it also indicates the *sign* of the result.
115 *
116 */
117 /*
118 * Constants:
119 * The hexadecimal values are the intended ones for the following
120 * constants. The decimal values may be used, provided that the
121 * compiler will convert from decimal to binary accurately enough
122 * to produce the hexadecimal values shown.
123 */
129 /*
130 * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
131 *
132 * integer array, contains the (24*i)-th to (24*i+23)-th
133 * bit of 2/pi after binary point. The corresponding
134 * floating value is
135 *
136 * ipio2[i] * 2^(-24(i+1)).
137 *
138 * NB: This table must have at least (e0-3)/24 + jk terms.
139 * For quad precision (e0 <= 16360, jk = 6), this is 686.
140 */
154 #if LDBL_MAX_EXP > 1024
259 #endif
260 };
271 };
274 {
278 /* initialize jk*/
282 /* determine jx,jv,q0, note that 3>q0 */
287 /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
292 /* compute q[0],q[1],...q[jk] */
297 }
300 recompute:
301 /* distill q[] into iq[] reversingly */
306 }
308 /* compute n */
318 }
330 }
333 }
340 }
341 }
346 }
347 }
349 /* check if recomputation is needed */
361 }
364 }
365 }
367 /* chop off zero terms */
372 jz--;
374 }
385 }
387 /* convert integer "bit" chunk to floating-point value */
392 }
394 /* compute PIo2[0,...,jp]*q[jz,...,0] */
399 }
401 /* compress fq[] into y[] */
414 // TODO: drop excess precision here once double_t is used
427 }
432 }
439 }
440 }
442 }