1 // SPDX-License-Identifier: GPL-2.0
2 /*---------------------------------------------------------------------------+
3 | poly_sin.c |
4 | |
5 | Computation of an approximation of the sin function and the cosine |
6 | function by a polynomial. |
7 | |
8 | Copyright (C) 1992,1993,1994,1997,1999 |
9 | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
10 | E-mail billm@melbpc.org.au |
11 | |
12 | |
13 +---------------------------------------------------------------------------*/
14
15 #include "exception.h"
16 #include "reg_constant.h"
17 #include "fpu_emu.h"
18 #include "fpu_system.h"
19 #include "control_w.h"
20 #include "poly.h"
21
22 #define N_COEFF_P 4
23 #define N_COEFF_N 4
24
25 static const unsigned long long pos_terms_l[N_COEFF_P] = {
26 0xaaaaaaaaaaaaaaabLL,
27 0x00d00d00d00cf906LL,
28 0x000006b99159a8bbLL,
29 0x000000000d7392e6LL
30 };
31
32 static const unsigned long long neg_terms_l[N_COEFF_N] = {
33 0x2222222222222167LL,
34 0x0002e3bc74aab624LL,
35 0x0000000b09229062LL,
36 0x00000000000c7973LL
37 };
38
39 #define N_COEFF_PH 4
40 #define N_COEFF_NH 4
41 static const unsigned long long pos_terms_h[N_COEFF_PH] = {
42 0x0000000000000000LL,
43 0x05b05b05b05b0406LL,
44 0x000049f93edd91a9LL,
45 0x00000000c9c9ed62LL
46 };
47
48 static const unsigned long long neg_terms_h[N_COEFF_NH] = {
49 0xaaaaaaaaaaaaaa98LL,
50 0x001a01a01a019064LL,
51 0x0000008f76c68a77LL,
52 0x0000000000d58f5eLL
53 };
54
55 /*--- poly_sine() -----------------------------------------------------------+
56 | |
57 +---------------------------------------------------------------------------*/
58 void poly_sine(FPU_REG *st0_ptr)
59 {
60 int exponent, echange;
61 Xsig accumulator, argSqrd, argTo4;
62 unsigned long fix_up, adj;
63 unsigned long long fixed_arg;
64 FPU_REG result;
65
66 exponent = exponent(st0_ptr);
67
68 accumulator.lsw = accumulator.midw = accumulator.msw = 0;
69
70 /* Split into two ranges, for arguments below and above 1.0 */
71 /* The boundary between upper and lower is approx 0.88309101259 */
72 if ((exponent < -1)
73 || ((exponent == -1) && (st0_ptr->sigh <= 0xe21240aa))) {
74 /* The argument is <= 0.88309101259 */
75
76 argSqrd.msw = st0_ptr->sigh;
77 argSqrd.midw = st0_ptr->sigl;
78 argSqrd.lsw = 0;
79 mul64_Xsig(&argSqrd, &significand(st0_ptr));
80 shr_Xsig(&argSqrd, 2 * (-1 - exponent));
81 argTo4.msw = argSqrd.msw;
82 argTo4.midw = argSqrd.midw;
83 argTo4.lsw = argSqrd.lsw;
84 mul_Xsig_Xsig(&argTo4, &argTo4);
85
86 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_l,
87 N_COEFF_N - 1);
88 mul_Xsig_Xsig(&accumulator, &argSqrd);
89 negate_Xsig(&accumulator);
90
91 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_l,
92 N_COEFF_P - 1);
93
94 shr_Xsig(&accumulator, 2); /* Divide by four */
95 accumulator.msw |= 0x80000000; /* Add 1.0 */
96
97 mul64_Xsig(&accumulator, &significand(st0_ptr));
98 mul64_Xsig(&accumulator, &significand(st0_ptr));
99 mul64_Xsig(&accumulator, &significand(st0_ptr));
100
101 /* Divide by four, FPU_REG compatible, etc */
102 exponent = 3 * exponent;
103
104 /* The minimum exponent difference is 3 */
105 shr_Xsig(&accumulator, exponent(st0_ptr) - exponent);
106
107 negate_Xsig(&accumulator);
108 XSIG_LL(accumulator) += significand(st0_ptr);
109
110 echange = round_Xsig(&accumulator);
111
112 setexponentpos(&result, exponent(st0_ptr) + echange);
113 } else {
114 /* The argument is > 0.88309101259 */
115 /* We use sin(st(0)) = cos(pi/2-st(0)) */
116
117 fixed_arg = significand(st0_ptr);
118
119 if (exponent == 0) {
120 /* The argument is >= 1.0 */
121
122 /* Put the binary point at the left. */
123 fixed_arg <<= 1;
124 }
125 /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
126 fixed_arg = 0x921fb54442d18469LL - fixed_arg;
127 /* There is a special case which arises due to rounding, to fix here. */
128 if (fixed_arg == 0xffffffffffffffffLL)
129 fixed_arg = 0;
130
131 XSIG_LL(argSqrd) = fixed_arg;
132 argSqrd.lsw = 0;
133 mul64_Xsig(&argSqrd, &fixed_arg);
134
135 XSIG_LL(argTo4) = XSIG_LL(argSqrd);
136 argTo4.lsw = argSqrd.lsw;
137 mul_Xsig_Xsig(&argTo4, &argTo4);
138
139 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_h,
140 N_COEFF_NH - 1);
141 mul_Xsig_Xsig(&accumulator, &argSqrd);
142 negate_Xsig(&accumulator);
143
144 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_h,
145 N_COEFF_PH - 1);
146 negate_Xsig(&accumulator);
147
148 mul64_Xsig(&accumulator, &fixed_arg);
149 mul64_Xsig(&accumulator, &fixed_arg);
150
151 shr_Xsig(&accumulator, 3);
152 negate_Xsig(&accumulator);
153
154 add_Xsig_Xsig(&accumulator, &argSqrd);
155
156 shr_Xsig(&accumulator, 1);
157
158 accumulator.lsw |= 1; /* A zero accumulator here would cause problems */
159 negate_Xsig(&accumulator);
160
161 /* The basic computation is complete. Now fix the answer to
162 compensate for the error due to the approximation used for
163 pi/2
164 */
165
166 /* This has an exponent of -65 */
167 fix_up = 0x898cc517;
168 /* The fix-up needs to be improved for larger args */
169 if (argSqrd.msw & 0xffc00000) {
170 /* Get about 32 bit precision in these: */
171 fix_up -= mul_32_32(0x898cc517, argSqrd.msw) / 6;
172 }
173 fix_up = mul_32_32(fix_up, LL_MSW(fixed_arg));
174
175 adj = accumulator.lsw; /* temp save */
176 accumulator.lsw -= fix_up;
177 if (accumulator.lsw > adj)
178 XSIG_LL(accumulator)--;
179
180 echange = round_Xsig(&accumulator);
181
182 setexponentpos(&result, echange - 1);
183 }
184
185 significand(&result) = XSIG_LL(accumulator);
186 setsign(&result, getsign(st0_ptr));
187 FPU_copy_to_reg0(&result, TAG_Valid);
188
189 #ifdef PARANOID
190 if ((exponent(&result) >= 0)
191 && (significand(&result) > 0x8000000000000000LL)) {
192 EXCEPTION(EX_INTERNAL | 0x150);
193 }
194 #endif /* PARANOID */
195
196 }
197
198 /*--- poly_cos() ------------------------------------------------------------+
199 | |
200 +---------------------------------------------------------------------------*/
201 void poly_cos(FPU_REG *st0_ptr)
202 {
203 FPU_REG result;
204 long int exponent, exp2, echange;
205 Xsig accumulator, argSqrd, fix_up, argTo4;
206 unsigned long long fixed_arg;
207
208 #ifdef PARANOID
209 if ((exponent(st0_ptr) > 0)
210 || ((exponent(st0_ptr) == 0)
211 && (significand(st0_ptr) > 0xc90fdaa22168c234LL))) {
212 EXCEPTION(EX_Invalid);
213 FPU_copy_to_reg0(&CONST_QNaN, TAG_Special);
214 return;
215 }
216 #endif /* PARANOID */
217
218 exponent = exponent(st0_ptr);
219
220 accumulator.lsw = accumulator.midw = accumulator.msw = 0;
221
222 if ((exponent < -1)
223 || ((exponent == -1) && (st0_ptr->sigh <= 0xb00d6f54))) {
224 /* arg is < 0.687705 */
225
226 argSqrd.msw = st0_ptr->sigh;
227 argSqrd.midw = st0_ptr->sigl;
228 argSqrd.lsw = 0;
229 mul64_Xsig(&argSqrd, &significand(st0_ptr));
230
231 if (exponent < -1) {
232 /* shift the argument right by the required places */
233 shr_Xsig(&argSqrd, 2 * (-1 - exponent));
234 }
235
236 argTo4.msw = argSqrd.msw;
237 argTo4.midw = argSqrd.midw;
238 argTo4.lsw = argSqrd.lsw;
239 mul_Xsig_Xsig(&argTo4, &argTo4);
240
241 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_h,
242 N_COEFF_NH - 1);
243 mul_Xsig_Xsig(&accumulator, &argSqrd);
244 negate_Xsig(&accumulator);
245
246 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_h,
247 N_COEFF_PH - 1);
248 negate_Xsig(&accumulator);
249
250 mul64_Xsig(&accumulator, &significand(st0_ptr));
251 mul64_Xsig(&accumulator, &significand(st0_ptr));
252 shr_Xsig(&accumulator, -2 * (1 + exponent));
253
254 shr_Xsig(&accumulator, 3);
255 negate_Xsig(&accumulator);
256
257 add_Xsig_Xsig(&accumulator, &argSqrd);
258
259 shr_Xsig(&accumulator, 1);
260
261 /* It doesn't matter if accumulator is all zero here, the
262 following code will work ok */
263 negate_Xsig(&accumulator);
264
265 if (accumulator.lsw & 0x80000000)
266 XSIG_LL(accumulator)++;
267 if (accumulator.msw == 0) {
268 /* The result is 1.0 */
269 FPU_copy_to_reg0(&CONST_1, TAG_Valid);
270 return;
271 } else {
272 significand(&result) = XSIG_LL(accumulator);
273
274 /* will be a valid positive nr with expon = -1 */
275 setexponentpos(&result, -1);
276 }
277 } else {
278 fixed_arg = significand(st0_ptr);
279
280 if (exponent == 0) {
281 /* The argument is >= 1.0 */
282
283 /* Put the binary point at the left. */
284 fixed_arg <<= 1;
285 }
286 /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
287 fixed_arg = 0x921fb54442d18469LL - fixed_arg;
288 /* There is a special case which arises due to rounding, to fix here. */
289 if (fixed_arg == 0xffffffffffffffffLL)
290 fixed_arg = 0;
291
292 exponent = -1;
293 exp2 = -1;
294
295 /* A shift is needed here only for a narrow range of arguments,
296 i.e. for fixed_arg approx 2^-32, but we pick up more... */
297 if (!(LL_MSW(fixed_arg) & 0xffff0000)) {
298 fixed_arg <<= 16;
299 exponent -= 16;
300 exp2 -= 16;
301 }
302
303 XSIG_LL(argSqrd) = fixed_arg;
304 argSqrd.lsw = 0;
305 mul64_Xsig(&argSqrd, &fixed_arg);
306
307 if (exponent < -1) {
308 /* shift the argument right by the required places */
309 shr_Xsig(&argSqrd, 2 * (-1 - exponent));
310 }
311
312 argTo4.msw = argSqrd.msw;
313 argTo4.midw = argSqrd.midw;
314 argTo4.lsw = argSqrd.lsw;
315 mul_Xsig_Xsig(&argTo4, &argTo4);
316
317 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_l,
318 N_COEFF_N - 1);
319 mul_Xsig_Xsig(&accumulator, &argSqrd);
320 negate_Xsig(&accumulator);
321
322 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_l,
323 N_COEFF_P - 1);
324
325 shr_Xsig(&accumulator, 2); /* Divide by four */
326 accumulator.msw |= 0x80000000; /* Add 1.0 */
327
328 mul64_Xsig(&accumulator, &fixed_arg);
329 mul64_Xsig(&accumulator, &fixed_arg);
330 mul64_Xsig(&accumulator, &fixed_arg);
331
332 /* Divide by four, FPU_REG compatible, etc */
333 exponent = 3 * exponent;
334
335 /* The minimum exponent difference is 3 */
336 shr_Xsig(&accumulator, exp2 - exponent);
337
338 negate_Xsig(&accumulator);
339 XSIG_LL(accumulator) += fixed_arg;
340
341 /* The basic computation is complete. Now fix the answer to
342 compensate for the error due to the approximation used for
343 pi/2
344 */
345
346 /* This has an exponent of -65 */
347 XSIG_LL(fix_up) = 0x898cc51701b839a2ll;
348 fix_up.lsw = 0;
349
350 /* The fix-up needs to be improved for larger args */
351 if (argSqrd.msw & 0xffc00000) {
352 /* Get about 32 bit precision in these: */
353 fix_up.msw -= mul_32_32(0x898cc517, argSqrd.msw) / 2;
354 fix_up.msw += mul_32_32(0x898cc517, argTo4.msw) / 24;
355 }
356
357 exp2 += norm_Xsig(&accumulator);
358 shr_Xsig(&accumulator, 1); /* Prevent overflow */
359 exp2++;
360 shr_Xsig(&fix_up, 65 + exp2);
361
362 add_Xsig_Xsig(&accumulator, &fix_up);
363
364 echange = round_Xsig(&accumulator);
365
366 setexponentpos(&result, exp2 + echange);
367 significand(&result) = XSIG_LL(accumulator);
368 }
369
370 FPU_copy_to_reg0(&result, TAG_Valid);
371
372 #ifdef PARANOID
373 if ((exponent(&result) >= 0)
374 && (significand(&result) > 0x8000000000000000LL)) {
375 EXCEPTION(EX_INTERNAL | 0x151);
376 }
377 #endif /* PARANOID */
378
379 }
380
Linux® is a registered trademark of Linus Torvalds in the United States and other countries.
TOMOYO® is a registered trademark of NTT DATA CORPORATION.