aboutsummaryrefslogtreecommitdiffstats
path: root/arch/x86/math-emu/poly_atan.c
diff options
context:
space:
mode:
Diffstat (limited to 'arch/x86/math-emu/poly_atan.c')
-rw-r--r--arch/x86/math-emu/poly_atan.c353
1 files changed, 166 insertions, 187 deletions
diff --git a/arch/x86/math-emu/poly_atan.c b/arch/x86/math-emu/poly_atan.c
index 82f702952f69..2f4ac8143fc3 100644
--- a/arch/x86/math-emu/poly_atan.c
+++ b/arch/x86/math-emu/poly_atan.c
@@ -18,28 +18,25 @@
18#include "control_w.h" 18#include "control_w.h"
19#include "poly.h" 19#include "poly.h"
20 20
21
22#define HIPOWERon 6 /* odd poly, negative terms */ 21#define HIPOWERon 6 /* odd poly, negative terms */
23static const unsigned long long oddnegterms[HIPOWERon] = 22static const unsigned long long oddnegterms[HIPOWERon] = {
24{ 23 0x0000000000000000LL, /* Dummy (not for - 1.0) */
25 0x0000000000000000LL, /* Dummy (not for - 1.0) */ 24 0x015328437f756467LL,
26 0x015328437f756467LL, 25 0x0005dda27b73dec6LL,
27 0x0005dda27b73dec6LL, 26 0x0000226bf2bfb91aLL,
28 0x0000226bf2bfb91aLL, 27 0x000000ccc439c5f7LL,
29 0x000000ccc439c5f7LL, 28 0x0000000355438407LL
30 0x0000000355438407LL 29};
31} ;
32 30
33#define HIPOWERop 6 /* odd poly, positive terms */ 31#define HIPOWERop 6 /* odd poly, positive terms */
34static const unsigned long long oddplterms[HIPOWERop] = 32static const unsigned long long oddplterms[HIPOWERop] = {
35{
36/* 0xaaaaaaaaaaaaaaabLL, transferred to fixedpterm[] */ 33/* 0xaaaaaaaaaaaaaaabLL, transferred to fixedpterm[] */
37 0x0db55a71875c9ac2LL, 34 0x0db55a71875c9ac2LL,
38 0x0029fce2d67880b0LL, 35 0x0029fce2d67880b0LL,
39 0x0000dfd3908b4596LL, 36 0x0000dfd3908b4596LL,
40 0x00000550fd61dab4LL, 37 0x00000550fd61dab4LL,
41 0x0000001c9422b3f9LL, 38 0x0000001c9422b3f9LL,
42 0x000000003e3301e1LL 39 0x000000003e3301e1LL
43}; 40};
44 41
45static const unsigned long long denomterm = 0xebd9b842c5c53a0eLL; 42static const unsigned long long denomterm = 0xebd9b842c5c53a0eLL;
@@ -48,182 +45,164 @@ static const Xsig fixedpterm = MK_XSIG(0xaaaaaaaa, 0xaaaaaaaa, 0xaaaaaaaa);
48 45
49static const Xsig pi_signif = MK_XSIG(0xc90fdaa2, 0x2168c234, 0xc4c6628b); 46static const Xsig pi_signif = MK_XSIG(0xc90fdaa2, 0x2168c234, 0xc4c6628b);
50 47
51
52/*--- poly_atan() -----------------------------------------------------------+ 48/*--- poly_atan() -----------------------------------------------------------+
53 | | 49 | |
54 +---------------------------------------------------------------------------*/ 50 +---------------------------------------------------------------------------*/
55void poly_atan(FPU_REG *st0_ptr, u_char st0_tag, 51void poly_atan(FPU_REG * st0_ptr, u_char st0_tag,
56 FPU_REG *st1_ptr, u_char st1_tag) 52 FPU_REG * st1_ptr, u_char st1_tag)
57{ 53{
58 u_char transformed, inverted, 54 u_char transformed, inverted, sign1, sign2;
59 sign1, sign2; 55 int exponent;
60 int exponent; 56 long int dummy_exp;
61 long int dummy_exp; 57 Xsig accumulator, Numer, Denom, accumulatore, argSignif, argSq, argSqSq;
62 Xsig accumulator, Numer, Denom, accumulatore, argSignif, 58 u_char tag;
63 argSq, argSqSq; 59
64 u_char tag; 60 sign1 = getsign(st0_ptr);
65 61 sign2 = getsign(st1_ptr);
66 sign1 = getsign(st0_ptr); 62 if (st0_tag == TAG_Valid) {
67 sign2 = getsign(st1_ptr); 63 exponent = exponent(st0_ptr);
68 if ( st0_tag == TAG_Valid ) 64 } else {
69 { 65 /* This gives non-compatible stack contents... */
70 exponent = exponent(st0_ptr); 66 FPU_to_exp16(st0_ptr, st0_ptr);
71 } 67 exponent = exponent16(st0_ptr);
72 else 68 }
73 { 69 if (st1_tag == TAG_Valid) {
74 /* This gives non-compatible stack contents... */ 70 exponent -= exponent(st1_ptr);
75 FPU_to_exp16(st0_ptr, st0_ptr); 71 } else {
76 exponent = exponent16(st0_ptr); 72 /* This gives non-compatible stack contents... */
77 } 73 FPU_to_exp16(st1_ptr, st1_ptr);
78 if ( st1_tag == TAG_Valid ) 74 exponent -= exponent16(st1_ptr);
79 { 75 }
80 exponent -= exponent(st1_ptr); 76
81 } 77 if ((exponent < 0) || ((exponent == 0) &&
82 else 78 ((st0_ptr->sigh < st1_ptr->sigh) ||
83 { 79 ((st0_ptr->sigh == st1_ptr->sigh) &&
84 /* This gives non-compatible stack contents... */ 80 (st0_ptr->sigl < st1_ptr->sigl))))) {
85 FPU_to_exp16(st1_ptr, st1_ptr); 81 inverted = 1;
86 exponent -= exponent16(st1_ptr); 82 Numer.lsw = Denom.lsw = 0;
87 } 83 XSIG_LL(Numer) = significand(st0_ptr);
88 84 XSIG_LL(Denom) = significand(st1_ptr);
89 if ( (exponent < 0) || ((exponent == 0) && 85 } else {
90 ((st0_ptr->sigh < st1_ptr->sigh) || 86 inverted = 0;
91 ((st0_ptr->sigh == st1_ptr->sigh) && 87 exponent = -exponent;
92 (st0_ptr->sigl < st1_ptr->sigl))) ) ) 88 Numer.lsw = Denom.lsw = 0;
93 { 89 XSIG_LL(Numer) = significand(st1_ptr);
94 inverted = 1; 90 XSIG_LL(Denom) = significand(st0_ptr);
95 Numer.lsw = Denom.lsw = 0; 91 }
96 XSIG_LL(Numer) = significand(st0_ptr); 92 div_Xsig(&Numer, &Denom, &argSignif);
97 XSIG_LL(Denom) = significand(st1_ptr); 93 exponent += norm_Xsig(&argSignif);
98 } 94
99 else 95 if ((exponent >= -1)
100 { 96 || ((exponent == -2) && (argSignif.msw > 0xd413ccd0))) {
101 inverted = 0; 97 /* The argument is greater than sqrt(2)-1 (=0.414213562...) */
102 exponent = -exponent; 98 /* Convert the argument by an identity for atan */
103 Numer.lsw = Denom.lsw = 0; 99 transformed = 1;
104 XSIG_LL(Numer) = significand(st1_ptr); 100
105 XSIG_LL(Denom) = significand(st0_ptr); 101 if (exponent >= 0) {
106 }
107 div_Xsig(&Numer, &Denom, &argSignif);
108 exponent += norm_Xsig(&argSignif);
109
110 if ( (exponent >= -1)
111 || ((exponent == -2) && (argSignif.msw > 0xd413ccd0)) )
112 {
113 /* The argument is greater than sqrt(2)-1 (=0.414213562...) */
114 /* Convert the argument by an identity for atan */
115 transformed = 1;
116
117 if ( exponent >= 0 )
118 {
119#ifdef PARANOID 102#ifdef PARANOID
120 if ( !( (exponent == 0) && 103 if (!((exponent == 0) &&
121 (argSignif.lsw == 0) && (argSignif.midw == 0) && 104 (argSignif.lsw == 0) && (argSignif.midw == 0) &&
122 (argSignif.msw == 0x80000000) ) ) 105 (argSignif.msw == 0x80000000))) {
123 { 106 EXCEPTION(EX_INTERNAL | 0x104); /* There must be a logic error */
124 EXCEPTION(EX_INTERNAL|0x104); /* There must be a logic error */ 107 return;
125 return; 108 }
126 }
127#endif /* PARANOID */ 109#endif /* PARANOID */
128 argSignif.msw = 0; /* Make the transformed arg -> 0.0 */ 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;
129 } 187 }
130 else 188
131 { 189 if (sign1) {
132 Numer.lsw = Denom.lsw = argSignif.lsw; 190 /* compute pi - accumulator */
133 XSIG_LL(Numer) = XSIG_LL(Denom) = XSIG_LL(argSignif); 191 shr_Xsig(&accumulator, 1 - exponent);
134 192 negate_Xsig(&accumulator);
135 if ( exponent < -1 ) 193 add_Xsig_Xsig(&accumulator, &pi_signif);
136 shr_Xsig(&Numer, -1-exponent); 194 exponent = 1;
137 negate_Xsig(&Numer);
138
139 shr_Xsig(&Denom, -exponent);
140 Denom.msw |= 0x80000000;
141
142 div_Xsig(&Numer, &Denom, &argSignif);
143
144 exponent = -1 + norm_Xsig(&argSignif);
145 } 195 }
146 } 196
147 else 197 exponent += round_Xsig(&accumulator);
148 { 198
149 transformed = 0; 199 significand(st1_ptr) = XSIG_LL(accumulator);
150 } 200 setexponent16(st1_ptr, exponent);
151 201
152 argSq.lsw = argSignif.lsw; argSq.midw = argSignif.midw; 202 tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign2);
153 argSq.msw = argSignif.msw; 203 FPU_settagi(1, tag);
154 mul_Xsig_Xsig(&argSq, &argSq); 204
155 205 set_precision_flag_up(); /* We do not really know if up or down,
156 argSqSq.lsw = argSq.lsw; argSqSq.midw = argSq.midw; argSqSq.msw = argSq.msw; 206 use this as the default. */
157 mul_Xsig_Xsig(&argSqSq, &argSqSq);
158
159 accumulatore.lsw = argSq.lsw;
160 XSIG_LL(accumulatore) = XSIG_LL(argSq);
161
162 shr_Xsig(&argSq, 2*(-1-exponent-1));
163 shr_Xsig(&argSqSq, 4*(-1-exponent-1));
164
165 /* Now have argSq etc with binary point at the left
166 .1xxxxxxxx */
167
168 /* Do the basic fixed point polynomial evaluation */
169 accumulator.msw = accumulator.midw = accumulator.lsw = 0;
170 polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq),
171 oddplterms, HIPOWERop-1);
172 mul64_Xsig(&accumulator, &XSIG_LL(argSq));
173 negate_Xsig(&accumulator);
174 polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq), oddnegterms, HIPOWERon-1);
175 negate_Xsig(&accumulator);
176 add_two_Xsig(&accumulator, &fixedpterm, &dummy_exp);
177
178 mul64_Xsig(&accumulatore, &denomterm);
179 shr_Xsig(&accumulatore, 1 + 2*(-1-exponent));
180 accumulatore.msw |= 0x80000000;
181
182 div_Xsig(&accumulator, &accumulatore, &accumulator);
183
184 mul_Xsig_Xsig(&accumulator, &argSignif);
185 mul_Xsig_Xsig(&accumulator, &argSq);
186
187 shr_Xsig(&accumulator, 3);
188 negate_Xsig(&accumulator);
189 add_Xsig_Xsig(&accumulator, &argSignif);
190
191 if ( transformed )
192 {
193 /* compute pi/4 - accumulator */
194 shr_Xsig(&accumulator, -1-exponent);
195 negate_Xsig(&accumulator);
196 add_Xsig_Xsig(&accumulator, &pi_signif);
197 exponent = -1;
198 }
199
200 if ( inverted )
201 {
202 /* compute pi/2 - accumulator */
203 shr_Xsig(&accumulator, -exponent);
204 negate_Xsig(&accumulator);
205 add_Xsig_Xsig(&accumulator, &pi_signif);
206 exponent = 0;
207 }
208
209 if ( sign1 )
210 {
211 /* compute pi - accumulator */
212 shr_Xsig(&accumulator, 1 - exponent);
213 negate_Xsig(&accumulator);
214 add_Xsig_Xsig(&accumulator, &pi_signif);
215 exponent = 1;
216 }
217
218 exponent += round_Xsig(&accumulator);
219
220 significand(st1_ptr) = XSIG_LL(accumulator);
221 setexponent16(st1_ptr, exponent);
222
223 tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign2);
224 FPU_settagi(1, tag);
225
226 set_precision_flag_up(); /* We do not really know if up or down,
227 use this as the default. */
228 207
229} 208}
generated by cgit v1.2.2 (git 2.25.0) at 2024-11-13 20:47:24 -0500