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