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