1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25
26 /*
27 * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31 /*
32 * __rem_pio2(x, y) passes back a better-than-double-precision
33 * approximation to x mod pi/2 in y[0]+y[1] and returns an integer
34 * congruent mod 8 to the integer part of x/(pi/2).
35 *
36 * This implementation tacitly assumes that x is finite and at
37 * least about pi/4 in magnitude.
38 */
39
40 #include "libm.h"
41
42 extern const int _TBL_ipio2_inf[];
43
44
45 /*
46 * invpio2: 53 bits of 2/pi
47 * pio2_1: first 33 bit of pi/2
48 * pio2_1t: pi/2 - pio2_1
49 * pio2_2: second 33 bit of pi/2
50 * pio2_2t: pi/2 - pio2_2
51 * pio2_3: third 33 bit of pi/2
52 * pio2_3t: pi/2 - pio2_3
53 */
54 static const double half = 0.5,
55 invpio2 = 0.636619772367581343075535, /* 2^ -1 * 1.45F306DC9C883 */
56 pio2_1 = 1.570796326734125614166, /* 2^ 0 * 1.921FB54400000 */
57 pio2_1t = 6.077100506506192601475e-11, /* 2^-34 * 1.0B4611A626331 */
58 pio2_2 = 6.077100506303965976596e-11, /* 2^-34 * 1.0B4611A600000 */
59 pio2_2t = 2.022266248795950732400e-21, /* 2^-69 * 1.3198A2E037073 */
60 pio2_3 = 2.022266248711166455796e-21, /* 2^-69 * 1.3198A2E000000 */
61 pio2_3t = 8.478427660368899643959e-32; /* 2^-104 * 1.B839A252049C1 */
62
63
64 int
65 __rem_pio2(double x, double *y)
66 {
67 double w, t, r, fn;
68 double tx[3];
69 int e0, i, j, nx, n, ix, hx, lx;
70
71 hx = ((int *)&x)[HIWORD];
72 ix = hx & 0x7fffffff;
73
74 if (ix < 0x4002d97c) {
75 /* |x| < 3pi/4, special case with n=1 */
76 t = fabs(x) - pio2_1;
77
78 if (ix != 0x3ff921fb) { /* 33+53 bit pi is good enough */
79 y[0] = t - pio2_1t;
80 y[1] = (t - y[0]) - pio2_1t;
81 } else { /* near pi/2, use 33+33+53 bit pi */
82 t -= pio2_2;
83 y[0] = t - pio2_2t;
84 y[1] = (t - y[0]) - pio2_2t;
85 }
86
87 if (hx < 0) {
88 y[0] = -y[0];
89 y[1] = -y[1];
90 return (-1);
91 }
92
93 return (1);
94 }
95
96 if (ix <= 0x413921fb) {
97 /* |x| <= 2^19 pi */
98 t = fabs(x);
99 n = (int)(t * invpio2 + half);
100 fn = (double)n;
101 r = t - fn * pio2_1;
102 j = ix >> 20;
103 w = fn * pio2_1t; /* 1st round good to 85 bit */
104 y[0] = r - w;
105 i = j - ((((int *)y)[HIWORD] >> 20) & 0x7ff);
106
107 if (i > 16) { /* 2nd iteration (rare) */
108 /* 2nd round good to 118 bit */
109 if (i < 35) {
110 t = r; /* r-fn*pio2_2 may not be exact */
111 w = fn * pio2_2;
112 r = t - w;
113 w = fn * pio2_2t - ((t - r) - w);
114 y[0] = r - w;
115 } else {
116 r -= fn * pio2_2;
117 w = fn * pio2_2t;
118 y[0] = r - w;
119 i = j - ((((int *)y)[HIWORD] >> 20) & 0x7ff);
120
121 if (i > 49) {
122 /* 3rd iteration (extremely rare) */
123 if (i < 68) {
124 t = r;
125 w = fn * pio2_3;
126 r = t - w;
127 w = fn * pio2_3t - ((t - r) -
128 w);
129 y[0] = r - w;
130 } else {
131 /*
132 * 3rd round good to 151 bits;
133 * covered all possible cases
134 */
135 r -= fn * pio2_3;
136 w = fn * pio2_3t;
137 y[0] = r - w;
138 }
139 }
140 }
141 }
142
143 y[1] = (r - y[0]) - w;
144
145 if (hx < 0) {
146 y[0] = -y[0];
147 y[1] = -y[1];
148 return (-n);
149 }
150
151 return (n);
152 }
153
154 e0 = (ix >> 20) - 1046; /* e0 = ilogb(x)-23; */
155
156 /* break x into three 24 bit pieces */
157 lx = ((int *)&x)[LOWORD];
158 i = (lx & 0x1f) << 19;
159 tx[2] = (double)i;
160 j = (lx >> 5) & 0xffffff;
161 tx[1] = (double)j;
162 tx[0] = (double)((((ix & 0xfffff) | 0x100000) << 3) | ((unsigned)lx >>
163 29));
164 nx = 3;
165
166 if (i == 0) {
167 /* skip zero term */
168 nx--;
169
170 if (j == 0)
171 nx--;
172 }
173
174 n = __rem_pio2m(tx, y, e0, nx, 2, _TBL_ipio2_inf);
175
176 if (hx < 0) {
177 y[0] = -y[0];
178 y[1] = -y[1];
179 return (-n);
180 }
181
182 return (n);
183 }