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
25static const unsigned long long pos_terms_l[N_COEFF_P] = {
26	0xaaaaaaaaaaaaaaabLL,
27	0x00d00d00d00cf906LL,
28	0x000006b99159a8bbLL,
29	0x000000000d7392e6LL
30};
31
32static 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
41static const unsigned long long pos_terms_h[N_COEFF_PH] = {
42	0x0000000000000000LL,
43	0x05b05b05b05b0406LL,
44	0x000049f93edd91a9LL,
45	0x00000000c9c9ed62LL
46};
47
48static 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 +---------------------------------------------------------------------------*/
58void 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 +---------------------------------------------------------------------------*/
201void 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