1 /* kf_rem_pio2.c -- float version of k_rem_pio2.c
2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
3 */
4
5 /*
6 * ====================================================
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8 *
9 * Developed at SunPro, a Sun Microsystems, Inc. business.
10 * Permission to use, copy, modify, and distribute this
11 * software is freely granted, provided that this notice
12 * is preserved.
13 * ====================================================
14 */
15
16 #include "fdlibm.h"
17
18 /* In the float version, the input parameter x contains 8 bit
19 integers, not 24 bit integers. 113 bit precision is not supported. */
20
21 static const int init_jk[] = { 4, 7, 9 }; /* initial value for jk */
22
23 static const float PIo2[] = {
24 1.5703125000e+00, /* 0x3fc90000 */
25 4.5776367188e-04, /* 0x39f00000 */
26 2.5987625122e-05, /* 0x37da0000 */
27 7.5437128544e-08, /* 0x33a20000 */
28 6.0026650317e-11, /* 0x2e840000 */
29 7.3896444519e-13, /* 0x2b500000 */
30 5.3845816694e-15, /* 0x27c20000 */
31 5.6378512969e-18, /* 0x22d00000 */
32 8.3009228831e-20, /* 0x1fc40000 */
33 3.2756352257e-22, /* 0x1bc60000 */
34 6.3331015649e-25, /* 0x17440000 */
35 };
36
37 static const float zero = 0.0, one = 1.0,
38 two8 = 2.5600000000e+02, /* 0x43800000 */
39 twon8 = 3.9062500000e-03; /* 0x3b800000 */
40
41 #ifdef __GNUC__
42 #pragma GCC diagnostic ignored "-Wpragmas"
43 #pragma GCC diagnostic ignored "-Wunknown-warning-option"
44 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
45 /* GCC analyzer gets confused about the use of 'iq' here */
46 #pragma GCC diagnostic ignored "-Wanalyzer-use-of-uninitialized-value"
47 #endif
48
49 int
__kernel_rem_pio2f(float * x,float * y,int e0,int nx,int prec,const __int32_t * ipio2)50 __kernel_rem_pio2f(float *x, float *y, int e0, int nx, int prec,
51 const __int32_t *ipio2)
52 {
53 __int32_t jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih;
54 float z, fw, f[20], fq[20], q[20];
55
56 /* initialize jk*/
57 jk = init_jk[prec];
58 jp = jk;
59
60 /* determine jx,jv,q0, note that 3>q0 */
61 jx = nx - 1;
62 jv = (e0 - 3) / 8;
63 if (jv < 0)
64 jv = 0;
65 q0 = e0 - 8 * (jv + 1);
66
67 /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
68 j = jv - jx;
69 m = jx + jk;
70 for (i = 0; i <= m; i++, j++)
71 f[i] = (j < 0) ? zero : (float)ipio2[j];
72
73 /* compute q[0],q[1],...q[jk] */
74 for (i = 0; i <= jk; i++) {
75 for (j = 0, fw = 0.0; j <= jx; j++)
76 fw += x[j] * f[jx + i - j];
77 q[i] = fw;
78 }
79
80 jz = jk;
81 recompute:
82 /* distill q[] into iq[] reversingly */
83 for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) {
84 fw = (float)((__int32_t)(twon8 * z));
85 iq[i] = (__int32_t)(z - two8 * fw);
86 z = q[j - 1] + fw;
87 }
88
89 /* compute n */
90 z = scalbnf(z, (int)q0); /* actual value of z */
91 z -= (float)8.0 * floorf(z * (float)0.125); /* trim off integer >= 8 */
92 n = (__int32_t)z;
93 z -= (float)n;
94 ih = 0;
95 if (q0 > 0) { /* need iq[jz-1] to determine n */
96 i = (iq[jz - 1] >> (8 - q0));
97 n += i;
98 iq[jz - 1] -= i << (8 - q0);
99 ih = iq[jz - 1] >> (7 - q0);
100 } else if (q0 == 0)
101 ih = iq[jz - 1] >> 7;
102 else if (z >= (float)0.5)
103 ih = 2;
104
105 if (ih > 0) { /* q > 0.5 */
106 n += 1;
107 carry = 0;
108 for (i = 0; i < jz; i++) { /* compute 1-q */
109 j = iq[i];
110 if (carry == 0) {
111 if (j != 0) {
112 carry = 1;
113 iq[i] = 0x100 - j;
114 }
115 } else
116 iq[i] = 0xff - j;
117 }
118 if (q0 > 0) { /* rare case: chance is 1 in 12 */
119 switch (q0) {
120 case 1:
121 iq[jz - 1] &= 0x7f;
122 break;
123 case 2:
124 iq[jz - 1] &= 0x3f;
125 break;
126 }
127 }
128 if (ih == 2) {
129 z = one - z;
130 if (carry != 0)
131 z -= scalbnf(one, (int)q0);
132 }
133 }
134
135 /* check if recomputation is needed */
136 if (z == zero) {
137 j = 0;
138 for (i = jz - 1; i >= jk; i--)
139 j |= iq[i];
140 if (j == 0) { /* need recomputation */
141 for (k = 1; iq[jk - k] == 0; k++)
142 ; /* k = no. of terms needed */
143
144 for (i = jz + 1; i <= jz + k; i++) { /* add q[jz+1] to q[jz+k] */
145 f[jx + i] = (float)ipio2[jv + i];
146 for (j = 0, fw = 0.0; j <= jx; j++)
147 fw += x[j] * f[jx + i - j];
148 q[i] = fw;
149 }
150 jz += k;
151 goto recompute;
152 }
153 }
154
155 /* chop off zero terms */
156 if (z == (float)0.0) {
157 jz -= 1;
158 q0 -= 8;
159 while (iq[jz] == 0) {
160 jz--;
161 q0 -= 8;
162 }
163 } else { /* break z into 8-bit if necessary */
164 z = scalbnf(z, -(int)q0);
165 if (z >= two8) {
166 fw = (float)((__int32_t)(twon8 * z));
167 iq[jz] = (__int32_t)(z - two8 * fw);
168 jz += 1;
169 q0 += 8;
170 iq[jz] = (__int32_t)fw;
171 } else
172 iq[jz] = (__int32_t)z;
173 }
174
175 /* convert integer "bit" chunk to floating-point value */
176 fw = scalbnf(one, (int)q0);
177 for (i = jz; i >= 0; i--) {
178 q[i] = fw * (float)iq[i];
179 fw *= twon8;
180 }
181
182 /* compute PIo2[0,...,jp]*q[jz,...,0] */
183 for (i = jz; i >= 0; i--) {
184 for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++)
185 fw += PIo2[k] * q[i + k];
186 fq[jz - i] = fw;
187 }
188
189 /* compress fq[] into y[] */
190 switch (prec) {
191 case 0:
192 fw = 0.0;
193 for (i = jz; i >= 0; i--)
194 fw += fq[i];
195 y[0] = (ih == 0) ? fw : -fw;
196 break;
197 case 1:
198 case 2:
199 fw = 0.0;
200 for (i = jz; i >= 0; i--)
201 fw += fq[i];
202 y[0] = (ih == 0) ? fw : -fw;
203 fw = fq[0] - fw;
204 for (i = 1; i <= jz; i++)
205 fw += fq[i];
206 y[1] = (ih == 0) ? fw : -fw;
207 break;
208 case 3: /* painful */
209 for (i = jz; i > 0; i--) {
210 fw = fq[i - 1] + fq[i];
211 fq[i] += fq[i - 1] - fw;
212 fq[i - 1] = fw;
213 }
214 for (i = jz; i > 1; i--) {
215 fw = fq[i - 1] + fq[i];
216 fq[i] += fq[i - 1] - fw;
217 fq[i - 1] = fw;
218 }
219 for (fw = 0.0, i = jz; i >= 2; i--)
220 fw += fq[i];
221 if (ih == 0) {
222 y[0] = fq[0];
223 y[1] = fq[1];
224 y[2] = fw;
225 } else {
226 y[0] = -fq[0];
227 y[1] = -fq[1];
228 y[2] = -fw;
229 }
230 }
231 return n & 7;
232 }
233
