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 #pragma GCC diagnostic ignored "-Wpragmas"
42 #pragma GCC diagnostic ignored "-Wunknown-warning-option"
43 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
44
45 int
__kernel_rem_pio2f(float * x,float * y,int e0,int nx,int prec,const __int32_t * ipio2)46 __kernel_rem_pio2f(float *x, float *y, int e0, int nx, int prec,
47 const __int32_t *ipio2)
48 {
49 __int32_t jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih;
50 float z, fw, f[20], fq[20], q[20];
51
52 /* initialize jk*/
53 jk = init_jk[prec];
54 jp = jk;
55
56 /* determine jx,jv,q0, note that 3>q0 */
57 jx = nx - 1;
58 jv = (e0 - 3) / 8;
59 if (jv < 0)
60 jv = 0;
61 q0 = e0 - 8 * (jv + 1);
62
63 /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
64 j = jv - jx;
65 m = jx + jk;
66 for (i = 0; i <= m; i++, j++)
67 f[i] = (j < 0) ? zero : (float)ipio2[j];
68
69 /* compute q[0],q[1],...q[jk] */
70 for (i = 0; i <= jk; i++) {
71 for (j = 0, fw = 0.0; j <= jx; j++)
72 fw += x[j] * f[jx + i - j];
73 q[i] = fw;
74 }
75
76 jz = jk;
77 recompute:
78 /* distill q[] into iq[] reversingly */
79 for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) {
80 fw = (float)((__int32_t)(twon8 * z));
81 iq[i] = (__int32_t)(z - two8 * fw);
82 z = q[j - 1] + fw;
83 }
84
85 /* compute n */
86 z = scalbnf(z, (int)q0); /* actual value of z */
87 z -= (float)8.0 * floorf(z * (float)0.125); /* trim off integer >= 8 */
88 n = (__int32_t)z;
89 z -= (float)n;
90 ih = 0;
91 if (q0 > 0) { /* need iq[jz-1] to determine n */
92 i = (iq[jz - 1] >> (8 - q0));
93 n += i;
94 iq[jz - 1] -= i << (8 - q0);
95 ih = iq[jz - 1] >> (7 - q0);
96 } else if (q0 == 0)
97 ih = iq[jz - 1] >> 7;
98 else if (z >= (float)0.5)
99 ih = 2;
100
101 if (ih > 0) { /* q > 0.5 */
102 n += 1;
103 carry = 0;
104 for (i = 0; i < jz; i++) { /* compute 1-q */
105 j = iq[i];
106 if (carry == 0) {
107 if (j != 0) {
108 carry = 1;
109 iq[i] = 0x100 - j;
110 }
111 } else
112 iq[i] = 0xff - j;
113 }
114 if (q0 > 0) { /* rare case: chance is 1 in 12 */
115 switch (q0) {
116 case 1:
117 iq[jz - 1] &= 0x7f;
118 break;
119 case 2:
120 iq[jz - 1] &= 0x3f;
121 break;
122 }
123 }
124 if (ih == 2) {
125 z = one - z;
126 if (carry != 0)
127 z -= scalbnf(one, (int)q0);
128 }
129 }
130
131 /* check if recomputation is needed */
132 if (z == zero) {
133 j = 0;
134 for (i = jz - 1; i >= jk; i--)
135 j |= iq[i];
136 if (j == 0) { /* need recomputation */
137 for (k = 1; iq[jk - k] == 0; k++)
138 ; /* k = no. of terms needed */
139
140 for (i = jz + 1; i <= jz + k; i++) { /* add q[jz+1] to q[jz+k] */
141 f[jx + i] = (float)ipio2[jv + i];
142 for (j = 0, fw = 0.0; j <= jx; j++)
143 fw += x[j] * f[jx + i - j];
144 q[i] = fw;
145 }
146 jz += k;
147 goto recompute;
148 }
149 }
150
151 /* chop off zero terms */
152 if (z == (float)0.0) {
153 jz -= 1;
154 q0 -= 8;
155 while (iq[jz] == 0) {
156 jz--;
157 q0 -= 8;
158 }
159 } else { /* break z into 8-bit if necessary */
160 z = scalbnf(z, -(int)q0);
161 if (z >= two8) {
162 fw = (float)((__int32_t)(twon8 * z));
163 iq[jz] = (__int32_t)(z - two8 * fw);
164 jz += 1;
165 q0 += 8;
166 iq[jz] = (__int32_t)fw;
167 } else
168 iq[jz] = (__int32_t)z;
169 }
170
171 /* convert integer "bit" chunk to floating-point value */
172 fw = scalbnf(one, (int)q0);
173 for (i = jz; i >= 0; i--) {
174 q[i] = fw * (float)iq[i];
175 fw *= twon8;
176 }
177
178 /* compute PIo2[0,...,jp]*q[jz,...,0] */
179 for (i = jz; i >= 0; i--) {
180 for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++)
181 fw += PIo2[k] * q[i + k];
182 fq[jz - i] = fw;
183 }
184
185 /* compress fq[] into y[] */
186 switch (prec) {
187 case 0:
188 fw = 0.0;
189 for (i = jz; i >= 0; i--)
190 fw += fq[i];
191 y[0] = (ih == 0) ? fw : -fw;
192 break;
193 case 1:
194 case 2:
195 fw = 0.0;
196 for (i = jz; i >= 0; i--)
197 fw += fq[i];
198 y[0] = (ih == 0) ? fw : -fw;
199 fw = fq[0] - fw;
200 for (i = 1; i <= jz; i++)
201 fw += fq[i];
202 y[1] = (ih == 0) ? fw : -fw;
203 break;
204 case 3: /* painful */
205 for (i = jz; i > 0; i--) {
206 fw = fq[i - 1] + fq[i];
207 fq[i] += fq[i - 1] - fw;
208 fq[i - 1] = fw;
209 }
210 for (i = jz; i > 1; i--) {
211 fw = fq[i - 1] + fq[i];
212 fq[i] += fq[i - 1] - fw;
213 fq[i - 1] = fw;
214 }
215 for (fw = 0.0, i = jz; i >= 2; i--)
216 fw += fq[i];
217 if (ih == 0) {
218 y[0] = fq[0];
219 y[1] = fq[1];
220 y[2] = fw;
221 } else {
222 y[0] = -fq[0];
223 y[1] = -fq[1];
224 y[2] = -fw;
225 }
226 }
227 return n & 7;
228 }
229
