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 */ |
23 | static 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 */ |
33 | static 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 | |
43 | static const unsigned long long denomterm = 0xebd9b842c5c53a0eLL; |
44 | |
45 | static const Xsig fixedpterm = MK_XSIG(0xaaaaaaaa, 0xaaaaaaaa, 0xaaaaaaaa); |
46 | |
47 | static const Xsig pi_signif = MK_XSIG(0xc90fdaa2, 0x2168c234, 0xc4c6628b); |
48 | |
49 | /*--- poly_atan() -----------------------------------------------------------+ |
50 | | | |
51 | +---------------------------------------------------------------------------*/ |
52 | void 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(a: st0_ptr, x: 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(a: st1_ptr, x: 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(x1: &Numer, x2: &Denom, dest: &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, n: -1 - exponent); |
118 | negate_Xsig(x: &Numer); |
119 | |
120 | shr_Xsig(&Denom, n: -exponent); |
121 | Denom.msw |= 0x80000000; |
122 | |
123 | div_Xsig(x1: &Numer, x2: &Denom, dest: &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(dest: &argSq, mult: &argSq); |
135 | |
136 | argSqSq.lsw = argSq.lsw; |
137 | argSqSq.midw = argSq.midw; |
138 | argSqSq.msw = argSq.msw; |
139 | mul_Xsig_Xsig(dest: &argSqSq, mult: &argSqSq); |
140 | |
141 | accumulatore.lsw = argSq.lsw; |
142 | XSIG_LL(accumulatore) = XSIG_LL(argSq); |
143 | |
144 | shr_Xsig(&argSq, n: 2 * (-1 - exponent - 1)); |
145 | shr_Xsig(&argSqSq, n: 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, x: &XSIG_LL(argSqSq), |
153 | terms: oddplterms, HIPOWERop - 1); |
154 | mul64_Xsig(&accumulator, mult: &XSIG_LL(argSq)); |
155 | negate_Xsig(x: &accumulator); |
156 | polynomial_Xsig(&accumulator, x: &XSIG_LL(argSqSq), terms: oddnegterms, |
157 | HIPOWERon - 1); |
158 | negate_Xsig(x: &accumulator); |
159 | add_two_Xsig(dest: &accumulator, x2: &fixedpterm, exp: &dummy_exp); |
160 | |
161 | mul64_Xsig(&accumulatore, mult: &denomterm); |
162 | shr_Xsig(&accumulatore, n: 1 + 2 * (-1 - exponent)); |
163 | accumulatore.msw |= 0x80000000; |
164 | |
165 | div_Xsig(x1: &accumulator, x2: &accumulatore, dest: &accumulator); |
166 | |
167 | mul_Xsig_Xsig(dest: &accumulator, mult: &argSignif); |
168 | mul_Xsig_Xsig(dest: &accumulator, mult: &argSq); |
169 | |
170 | shr_Xsig(&accumulator, n: 3); |
171 | negate_Xsig(x: &accumulator); |
172 | add_Xsig_Xsig(dest: &accumulator, x2: &argSignif); |
173 | |
174 | if (transformed) { |
175 | /* compute pi/4 - accumulator */ |
176 | shr_Xsig(&accumulator, n: -1 - exponent); |
177 | negate_Xsig(x: &accumulator); |
178 | add_Xsig_Xsig(dest: &accumulator, x2: &pi_signif); |
179 | exponent = -1; |
180 | } |
181 | |
182 | if (inverted) { |
183 | /* compute pi/2 - accumulator */ |
184 | shr_Xsig(&accumulator, n: -exponent); |
185 | negate_Xsig(x: &accumulator); |
186 | add_Xsig_Xsig(dest: &accumulator, x2: &pi_signif); |
187 | exponent = 0; |
188 | } |
189 | |
190 | if (sign1) { |
191 | /* compute pi - accumulator */ |
192 | shr_Xsig(&accumulator, n: 1 - exponent); |
193 | negate_Xsig(x: &accumulator); |
194 | add_Xsig_Xsig(dest: &accumulator, x2: &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(arg: st1_ptr, extent: 1, dummy: 0, FULL_PRECISION, sign: sign2); |
204 | FPU_settagi(stnr: 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 | |