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