1// SPDX-License-Identifier: GPL-2.0
2/*---------------------------------------------------------------------------+
3 |  poly_atan.c                                                              |
4 |                                                                           |
5 | Compute the arctan of a FPU_REG, using a polynomial approximation.        |
6 |                                                                           |
7 | Copyright (C) 1992,1993,1994,1997                                         |
8 |                  W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
9 |                  E-mail   billm@suburbia.net                              |
10 |                                                                           |
11 |                                                                           |
12 +---------------------------------------------------------------------------*/
13
14#include "exception.h"
15#include "reg_constant.h"
16#include "fpu_emu.h"
17#include "fpu_system.h"
18#include "status_w.h"
19#include "control_w.h"
20#include "poly.h"
21
22#define	HIPOWERon	6	/* odd poly, negative terms */
23static const unsigned long long oddnegterms[HIPOWERon] = {
24	0x0000000000000000LL,	/* Dummy (not for - 1.0) */
25	0x015328437f756467LL,
26	0x0005dda27b73dec6LL,
27	0x0000226bf2bfb91aLL,
28	0x000000ccc439c5f7LL,
29	0x0000000355438407LL
30};
31
32#define	HIPOWERop	6	/* odd poly, positive terms */
33static const unsigned long long oddplterms[HIPOWERop] = {
34/*  0xaaaaaaaaaaaaaaabLL,  transferred to fixedpterm[] */
35	0x0db55a71875c9ac2LL,
36	0x0029fce2d67880b0LL,
37	0x0000dfd3908b4596LL,
38	0x00000550fd61dab4LL,
39	0x0000001c9422b3f9LL,
40	0x000000003e3301e1LL
41};
42
43static const unsigned long long denomterm = 0xebd9b842c5c53a0eLL;
44
45static const Xsig fixedpterm = MK_XSIG(0xaaaaaaaa, 0xaaaaaaaa, 0xaaaaaaaa);
46
47static const Xsig pi_signif = MK_XSIG(0xc90fdaa2, 0x2168c234, 0xc4c6628b);
48
49/*--- poly_atan() -----------------------------------------------------------+
50 |                                                                           |
51 +---------------------------------------------------------------------------*/
52void poly_atan(FPU_REG *st0_ptr, u_char st0_tag,
53	       FPU_REG *st1_ptr, u_char st1_tag)
54{
55	u_char transformed, inverted, sign1, sign2;
56	int exponent;
57	long int dummy_exp;
58	Xsig accumulator, Numer, Denom, accumulatore, argSignif, argSq, argSqSq;
59	u_char tag;
60
61	sign1 = getsign(st0_ptr);
62	sign2 = getsign(st1_ptr);
63	if (st0_tag == TAG_Valid) {
64		exponent = exponent(st0_ptr);
65	} else {
66		/* This gives non-compatible stack contents... */
67		FPU_to_exp16(st0_ptr, st0_ptr);
68		exponent = exponent16(st0_ptr);
69	}
70	if (st1_tag == TAG_Valid) {
71		exponent -= exponent(st1_ptr);
72	} else {
73		/* This gives non-compatible stack contents... */
74		FPU_to_exp16(st1_ptr, st1_ptr);
75		exponent -= exponent16(st1_ptr);
76	}
77
78	if ((exponent < 0) || ((exponent == 0) &&
79			       ((st0_ptr->sigh < st1_ptr->sigh) ||
80				((st0_ptr->sigh == st1_ptr->sigh) &&
81				 (st0_ptr->sigl < st1_ptr->sigl))))) {
82		inverted = 1;
83		Numer.lsw = Denom.lsw = 0;
84		XSIG_LL(Numer) = significand(st0_ptr);
85		XSIG_LL(Denom) = significand(st1_ptr);
86	} else {
87		inverted = 0;
88		exponent = -exponent;
89		Numer.lsw = Denom.lsw = 0;
90		XSIG_LL(Numer) = significand(st1_ptr);
91		XSIG_LL(Denom) = significand(st0_ptr);
92	}
93	div_Xsig(&Numer, &Denom, &argSignif);
94	exponent += norm_Xsig(&argSignif);
95
96	if ((exponent >= -1)
97	    || ((exponent == -2) && (argSignif.msw > 0xd413ccd0))) {
98		/* The argument is greater than sqrt(2)-1 (=0.414213562...) */
99		/* Convert the argument by an identity for atan */
100		transformed = 1;
101
102		if (exponent >= 0) {
103#ifdef PARANOID
104			if (!((exponent == 0) &&
105			      (argSignif.lsw == 0) && (argSignif.midw == 0) &&
106			      (argSignif.msw == 0x80000000))) {
107				EXCEPTION(EX_INTERNAL | 0x104);	/* There must be a logic error */
108				return;
109			}
110#endif /* PARANOID */
111			argSignif.msw = 0;	/* Make the transformed arg -> 0.0 */
112		} else {
113			Numer.lsw = Denom.lsw = argSignif.lsw;
114			XSIG_LL(Numer) = XSIG_LL(Denom) = XSIG_LL(argSignif);
115
116			if (exponent < -1)
117				shr_Xsig(&Numer, -1 - exponent);
118			negate_Xsig(&Numer);
119
120			shr_Xsig(&Denom, -exponent);
121			Denom.msw |= 0x80000000;
122
123			div_Xsig(&Numer, &Denom, &argSignif);
124
125			exponent = -1 + norm_Xsig(&argSignif);
126		}
127	} else {
128		transformed = 0;
129	}
130
131	argSq.lsw = argSignif.lsw;
132	argSq.midw = argSignif.midw;
133	argSq.msw = argSignif.msw;
134	mul_Xsig_Xsig(&argSq, &argSq);
135
136	argSqSq.lsw = argSq.lsw;
137	argSqSq.midw = argSq.midw;
138	argSqSq.msw = argSq.msw;
139	mul_Xsig_Xsig(&argSqSq, &argSqSq);
140
141	accumulatore.lsw = argSq.lsw;
142	XSIG_LL(accumulatore) = XSIG_LL(argSq);
143
144	shr_Xsig(&argSq, 2 * (-1 - exponent - 1));
145	shr_Xsig(&argSqSq, 4 * (-1 - exponent - 1));
146
147	/* Now have argSq etc with binary point at the left
148	   .1xxxxxxxx */
149
150	/* Do the basic fixed point polynomial evaluation */
151	accumulator.msw = accumulator.midw = accumulator.lsw = 0;
152	polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq),
153			oddplterms, HIPOWERop - 1);
154	mul64_Xsig(&accumulator, &XSIG_LL(argSq));
155	negate_Xsig(&accumulator);
156	polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq), oddnegterms,
157			HIPOWERon - 1);
158	negate_Xsig(&accumulator);
159	add_two_Xsig(&accumulator, &fixedpterm, &dummy_exp);
160
161	mul64_Xsig(&accumulatore, &denomterm);
162	shr_Xsig(&accumulatore, 1 + 2 * (-1 - exponent));
163	accumulatore.msw |= 0x80000000;
164
165	div_Xsig(&accumulator, &accumulatore, &accumulator);
166
167	mul_Xsig_Xsig(&accumulator, &argSignif);
168	mul_Xsig_Xsig(&accumulator, &argSq);
169
170	shr_Xsig(&accumulator, 3);
171	negate_Xsig(&accumulator);
172	add_Xsig_Xsig(&accumulator, &argSignif);
173
174	if (transformed) {
175		/* compute pi/4 - accumulator */
176		shr_Xsig(&accumulator, -1 - exponent);
177		negate_Xsig(&accumulator);
178		add_Xsig_Xsig(&accumulator, &pi_signif);
179		exponent = -1;
180	}
181
182	if (inverted) {
183		/* compute pi/2 - accumulator */
184		shr_Xsig(&accumulator, -exponent);
185		negate_Xsig(&accumulator);
186		add_Xsig_Xsig(&accumulator, &pi_signif);
187		exponent = 0;
188	}
189
190	if (sign1) {
191		/* compute pi - accumulator */
192		shr_Xsig(&accumulator, 1 - exponent);
193		negate_Xsig(&accumulator);
194		add_Xsig_Xsig(&accumulator, &pi_signif);
195		exponent = 1;
196	}
197
198	exponent += round_Xsig(&accumulator);
199
200	significand(st1_ptr) = XSIG_LL(accumulator);
201	setexponent16(st1_ptr, exponent);
202
203	tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign2);
204	FPU_settagi(1, tag);
205
206	set_precision_flag_up();	/* We do not really know if up or down,
207					   use this as the default. */
208
209}
210