aboutsummaryrefslogtreecommitdiffstats
path: root/arch/x86/math-emu/poly_l2.c
diff options
context:
space:
mode:
Diffstat (limited to 'arch/x86/math-emu/poly_l2.c')
-rw-r--r--arch/x86/math-emu/poly_l2.c378
1 files changed, 175 insertions, 203 deletions
diff --git a/arch/x86/math-emu/poly_l2.c b/arch/x86/math-emu/poly_l2.c
index dd00e1d5b074..c0102ae87511 100644
--- a/arch/x86/math-emu/poly_l2.c
+++ b/arch/x86/math-emu/poly_l2.c
@@ -10,7 +10,6 @@
10 | | 10 | |
11 +---------------------------------------------------------------------------*/ 11 +---------------------------------------------------------------------------*/
12 12
13
14#include "exception.h" 13#include "exception.h"
15#include "reg_constant.h" 14#include "reg_constant.h"
16#include "fpu_emu.h" 15#include "fpu_emu.h"
@@ -18,255 +17,228 @@
18#include "control_w.h" 17#include "control_w.h"
19#include "poly.h" 18#include "poly.h"
20 19
21
22static void log2_kernel(FPU_REG const *arg, u_char argsign, 20static void log2_kernel(FPU_REG const *arg, u_char argsign,
23 Xsig *accum_result, long int *expon); 21 Xsig * accum_result, long int *expon);
24
25 22
26/*--- poly_l2() -------------------------------------------------------------+ 23/*--- poly_l2() -------------------------------------------------------------+
27 | Base 2 logarithm by a polynomial approximation. | 24 | Base 2 logarithm by a polynomial approximation. |
28 +---------------------------------------------------------------------------*/ 25 +---------------------------------------------------------------------------*/
29void poly_l2(FPU_REG *st0_ptr, FPU_REG *st1_ptr, u_char st1_sign) 26void poly_l2(FPU_REG * st0_ptr, FPU_REG * st1_ptr, u_char st1_sign)
30{ 27{
31 long int exponent, expon, expon_expon; 28 long int exponent, expon, expon_expon;
32 Xsig accumulator, expon_accum, yaccum; 29 Xsig accumulator, expon_accum, yaccum;
33 u_char sign, argsign; 30 u_char sign, argsign;
34 FPU_REG x; 31 FPU_REG x;
35 int tag; 32 int tag;
36 33
37 exponent = exponent16(st0_ptr); 34 exponent = exponent16(st0_ptr);
38 35
39 /* From st0_ptr, make a number > sqrt(2)/2 and < sqrt(2) */ 36 /* From st0_ptr, make a number > sqrt(2)/2 and < sqrt(2) */
40 if ( st0_ptr->sigh > (unsigned)0xb504f334 ) 37 if (st0_ptr->sigh > (unsigned)0xb504f334) {
41 { 38 /* Treat as sqrt(2)/2 < st0_ptr < 1 */
42 /* Treat as sqrt(2)/2 < st0_ptr < 1 */ 39 significand(&x) = -significand(st0_ptr);
43 significand(&x) = - significand(st0_ptr); 40 setexponent16(&x, -1);
44 setexponent16(&x, -1); 41 exponent++;
45 exponent++; 42 argsign = SIGN_NEG;
46 argsign = SIGN_NEG; 43 } else {
47 } 44 /* Treat as 1 <= st0_ptr < sqrt(2) */
48 else 45 x.sigh = st0_ptr->sigh - 0x80000000;
49 { 46 x.sigl = st0_ptr->sigl;
50 /* Treat as 1 <= st0_ptr < sqrt(2) */ 47 setexponent16(&x, 0);
51 x.sigh = st0_ptr->sigh - 0x80000000; 48 argsign = SIGN_POS;
52 x.sigl = st0_ptr->sigl; 49 }
53 setexponent16(&x, 0); 50 tag = FPU_normalize_nuo(&x);
54 argsign = SIGN_POS;
55 }
56 tag = FPU_normalize_nuo(&x);
57
58 if ( tag == TAG_Zero )
59 {
60 expon = 0;
61 accumulator.msw = accumulator.midw = accumulator.lsw = 0;
62 }
63 else
64 {
65 log2_kernel(&x, argsign, &accumulator, &expon);
66 }
67
68 if ( exponent < 0 )
69 {
70 sign = SIGN_NEG;
71 exponent = -exponent;
72 }
73 else
74 sign = SIGN_POS;
75 expon_accum.msw = exponent; expon_accum.midw = expon_accum.lsw = 0;
76 if ( exponent )
77 {
78 expon_expon = 31 + norm_Xsig(&expon_accum);
79 shr_Xsig(&accumulator, expon_expon - expon);
80
81 if ( sign ^ argsign )
82 negate_Xsig(&accumulator);
83 add_Xsig_Xsig(&accumulator, &expon_accum);
84 }
85 else
86 {
87 expon_expon = expon;
88 sign = argsign;
89 }
90
91 yaccum.lsw = 0; XSIG_LL(yaccum) = significand(st1_ptr);
92 mul_Xsig_Xsig(&accumulator, &yaccum);
93
94 expon_expon += round_Xsig(&accumulator);
95
96 if ( accumulator.msw == 0 )
97 {
98 FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
99 return;
100 }
101
102 significand(st1_ptr) = XSIG_LL(accumulator);
103 setexponent16(st1_ptr, expon_expon + exponent16(st1_ptr) + 1);
104
105 tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign ^ st1_sign);
106 FPU_settagi(1, tag);
107
108 set_precision_flag_up(); /* 80486 appears to always do this */
109
110 return;
111 51
112} 52 if (tag == TAG_Zero) {
53 expon = 0;
54 accumulator.msw = accumulator.midw = accumulator.lsw = 0;
55 } else {
56 log2_kernel(&x, argsign, &accumulator, &expon);
57 }
58
59 if (exponent < 0) {
60 sign = SIGN_NEG;
61 exponent = -exponent;
62 } else
63 sign = SIGN_POS;
64 expon_accum.msw = exponent;
65 expon_accum.midw = expon_accum.lsw = 0;
66 if (exponent) {
67 expon_expon = 31 + norm_Xsig(&expon_accum);
68 shr_Xsig(&accumulator, expon_expon - expon);
69
70 if (sign ^ argsign)
71 negate_Xsig(&accumulator);
72 add_Xsig_Xsig(&accumulator, &expon_accum);
73 } else {
74 expon_expon = expon;
75 sign = argsign;
76 }
77
78 yaccum.lsw = 0;
79 XSIG_LL(yaccum) = significand(st1_ptr);
80 mul_Xsig_Xsig(&accumulator, &yaccum);
81
82 expon_expon += round_Xsig(&accumulator);
83
84 if (accumulator.msw == 0) {
85 FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
86 return;
87 }
88
89 significand(st1_ptr) = XSIG_LL(accumulator);
90 setexponent16(st1_ptr, expon_expon + exponent16(st1_ptr) + 1);
113 91
92 tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign ^ st1_sign);
93 FPU_settagi(1, tag);
94
95 set_precision_flag_up(); /* 80486 appears to always do this */
96
97 return;
98
99}
114 100
115/*--- poly_l2p1() -----------------------------------------------------------+ 101/*--- poly_l2p1() -----------------------------------------------------------+
116 | Base 2 logarithm by a polynomial approximation. | 102 | Base 2 logarithm by a polynomial approximation. |
117 | log2(x+1) | 103 | log2(x+1) |
118 +---------------------------------------------------------------------------*/ 104 +---------------------------------------------------------------------------*/
119int poly_l2p1(u_char sign0, u_char sign1, 105int poly_l2p1(u_char sign0, u_char sign1,
120 FPU_REG *st0_ptr, FPU_REG *st1_ptr, FPU_REG *dest) 106 FPU_REG * st0_ptr, FPU_REG * st1_ptr, FPU_REG * dest)
121{ 107{
122 u_char tag; 108 u_char tag;
123 long int exponent; 109 long int exponent;
124 Xsig accumulator, yaccum; 110 Xsig accumulator, yaccum;
125 111
126 if ( exponent16(st0_ptr) < 0 ) 112 if (exponent16(st0_ptr) < 0) {
127 { 113 log2_kernel(st0_ptr, sign0, &accumulator, &exponent);
128 log2_kernel(st0_ptr, sign0, &accumulator, &exponent);
129 114
130 yaccum.lsw = 0; 115 yaccum.lsw = 0;
131 XSIG_LL(yaccum) = significand(st1_ptr); 116 XSIG_LL(yaccum) = significand(st1_ptr);
132 mul_Xsig_Xsig(&accumulator, &yaccum); 117 mul_Xsig_Xsig(&accumulator, &yaccum);
133 118
134 exponent += round_Xsig(&accumulator); 119 exponent += round_Xsig(&accumulator);
135 120
136 exponent += exponent16(st1_ptr) + 1; 121 exponent += exponent16(st1_ptr) + 1;
137 if ( exponent < EXP_WAY_UNDER ) exponent = EXP_WAY_UNDER; 122 if (exponent < EXP_WAY_UNDER)
123 exponent = EXP_WAY_UNDER;
138 124
139 significand(dest) = XSIG_LL(accumulator); 125 significand(dest) = XSIG_LL(accumulator);
140 setexponent16(dest, exponent); 126 setexponent16(dest, exponent);
141 127
142 tag = FPU_round(dest, 1, 0, FULL_PRECISION, sign0 ^ sign1); 128 tag = FPU_round(dest, 1, 0, FULL_PRECISION, sign0 ^ sign1);
143 FPU_settagi(1, tag); 129 FPU_settagi(1, tag);
144 130
145 if ( tag == TAG_Valid ) 131 if (tag == TAG_Valid)
146 set_precision_flag_up(); /* 80486 appears to always do this */ 132 set_precision_flag_up(); /* 80486 appears to always do this */
147 } 133 } else {
148 else 134 /* The magnitude of st0_ptr is far too large. */
149 {
150 /* The magnitude of st0_ptr is far too large. */
151 135
152 if ( sign0 != SIGN_POS ) 136 if (sign0 != SIGN_POS) {
153 { 137 /* Trying to get the log of a negative number. */
154 /* Trying to get the log of a negative number. */ 138#ifdef PECULIAR_486 /* Stupid 80486 doesn't worry about log(negative). */
155#ifdef PECULIAR_486 /* Stupid 80486 doesn't worry about log(negative). */ 139 changesign(st1_ptr);
156 changesign(st1_ptr);
157#else 140#else
158 if ( arith_invalid(1) < 0 ) 141 if (arith_invalid(1) < 0)
159 return 1; 142 return 1;
160#endif /* PECULIAR_486 */ 143#endif /* PECULIAR_486 */
161 } 144 }
162 145
163 /* 80486 appears to do this */ 146 /* 80486 appears to do this */
164 if ( sign0 == SIGN_NEG ) 147 if (sign0 == SIGN_NEG)
165 set_precision_flag_down(); 148 set_precision_flag_down();
166 else 149 else
167 set_precision_flag_up(); 150 set_precision_flag_up();
168 } 151 }
169 152
170 if ( exponent(dest) <= EXP_UNDER ) 153 if (exponent(dest) <= EXP_UNDER)
171 EXCEPTION(EX_Underflow); 154 EXCEPTION(EX_Underflow);
172 155
173 return 0; 156 return 0;
174 157
175} 158}
176 159
177
178
179
180#undef HIPOWER 160#undef HIPOWER
181#define HIPOWER 10 161#define HIPOWER 10
182static const unsigned long long logterms[HIPOWER] = 162static const unsigned long long logterms[HIPOWER] = {
183{ 163 0x2a8eca5705fc2ef0LL,
184 0x2a8eca5705fc2ef0LL, 164 0xf6384ee1d01febceLL,
185 0xf6384ee1d01febceLL, 165 0x093bb62877cdf642LL,
186 0x093bb62877cdf642LL, 166 0x006985d8a9ec439bLL,
187 0x006985d8a9ec439bLL, 167 0x0005212c4f55a9c8LL,
188 0x0005212c4f55a9c8LL, 168 0x00004326a16927f0LL,
189 0x00004326a16927f0LL, 169 0x0000038d1d80a0e7LL,
190 0x0000038d1d80a0e7LL, 170 0x0000003141cc80c6LL,
191 0x0000003141cc80c6LL, 171 0x00000002b1668c9fLL,
192 0x00000002b1668c9fLL, 172 0x000000002c7a46aaLL
193 0x000000002c7a46aaLL
194}; 173};
195 174
196static const unsigned long leadterm = 0xb8000000; 175static const unsigned long leadterm = 0xb8000000;
197 176
198
199/*--- log2_kernel() ---------------------------------------------------------+ 177/*--- log2_kernel() ---------------------------------------------------------+
200 | Base 2 logarithm by a polynomial approximation. | 178 | Base 2 logarithm by a polynomial approximation. |
201 | log2(x+1) | 179 | log2(x+1) |
202 +---------------------------------------------------------------------------*/ 180 +---------------------------------------------------------------------------*/
203static void log2_kernel(FPU_REG const *arg, u_char argsign, Xsig *accum_result, 181static void log2_kernel(FPU_REG const *arg, u_char argsign, Xsig * accum_result,
204 long int *expon) 182 long int *expon)
205{ 183{
206 long int exponent, adj; 184 long int exponent, adj;
207 unsigned long long Xsq; 185 unsigned long long Xsq;
208 Xsig accumulator, Numer, Denom, argSignif, arg_signif; 186 Xsig accumulator, Numer, Denom, argSignif, arg_signif;
209 187
210 exponent = exponent16(arg); 188 exponent = exponent16(arg);
211 Numer.lsw = Denom.lsw = 0; 189 Numer.lsw = Denom.lsw = 0;
212 XSIG_LL(Numer) = XSIG_LL(Denom) = significand(arg); 190 XSIG_LL(Numer) = XSIG_LL(Denom) = significand(arg);
213 if ( argsign == SIGN_POS ) 191 if (argsign == SIGN_POS) {
214 { 192 shr_Xsig(&Denom, 2 - (1 + exponent));
215 shr_Xsig(&Denom, 2 - (1 + exponent)); 193 Denom.msw |= 0x80000000;
216 Denom.msw |= 0x80000000; 194 div_Xsig(&Numer, &Denom, &argSignif);
217 div_Xsig(&Numer, &Denom, &argSignif); 195 } else {
218 } 196 shr_Xsig(&Denom, 1 - (1 + exponent));
219 else 197 negate_Xsig(&Denom);
220 { 198 if (Denom.msw & 0x80000000) {
221 shr_Xsig(&Denom, 1 - (1 + exponent)); 199 div_Xsig(&Numer, &Denom, &argSignif);
222 negate_Xsig(&Denom); 200 exponent++;
223 if ( Denom.msw & 0x80000000 ) 201 } else {
224 { 202 /* Denom must be 1.0 */
225 div_Xsig(&Numer, &Denom, &argSignif); 203 argSignif.lsw = Numer.lsw;
226 exponent ++; 204 argSignif.midw = Numer.midw;
227 } 205 argSignif.msw = Numer.msw;
228 else 206 }
229 {
230 /* Denom must be 1.0 */
231 argSignif.lsw = Numer.lsw; argSignif.midw = Numer.midw;
232 argSignif.msw = Numer.msw;
233 } 207 }
234 }
235 208
236#ifndef PECULIAR_486 209#ifndef PECULIAR_486
237 /* Should check here that |local_arg| is within the valid range */ 210 /* Should check here that |local_arg| is within the valid range */
238 if ( exponent >= -2 ) 211 if (exponent >= -2) {
239 { 212 if ((exponent > -2) || (argSignif.msw > (unsigned)0xafb0ccc0)) {
240 if ( (exponent > -2) || 213 /* The argument is too large */
241 (argSignif.msw > (unsigned)0xafb0ccc0) ) 214 }
242 {
243 /* The argument is too large */
244 } 215 }
245 }
246#endif /* PECULIAR_486 */ 216#endif /* PECULIAR_486 */
247 217
248 arg_signif.lsw = argSignif.lsw; XSIG_LL(arg_signif) = XSIG_LL(argSignif); 218 arg_signif.lsw = argSignif.lsw;
249 adj = norm_Xsig(&argSignif); 219 XSIG_LL(arg_signif) = XSIG_LL(argSignif);
250 accumulator.lsw = argSignif.lsw; XSIG_LL(accumulator) = XSIG_LL(argSignif); 220 adj = norm_Xsig(&argSignif);
251 mul_Xsig_Xsig(&accumulator, &accumulator); 221 accumulator.lsw = argSignif.lsw;
252 shr_Xsig(&accumulator, 2*(-1 - (1 + exponent + adj))); 222 XSIG_LL(accumulator) = XSIG_LL(argSignif);
253 Xsq = XSIG_LL(accumulator); 223 mul_Xsig_Xsig(&accumulator, &accumulator);
254 if ( accumulator.lsw & 0x80000000 ) 224 shr_Xsig(&accumulator, 2 * (-1 - (1 + exponent + adj)));
255 Xsq++; 225 Xsq = XSIG_LL(accumulator);
256 226 if (accumulator.lsw & 0x80000000)
257 accumulator.msw = accumulator.midw = accumulator.lsw = 0; 227 Xsq++;
258 /* Do the basic fixed point polynomial evaluation */ 228
259 polynomial_Xsig(&accumulator, &Xsq, logterms, HIPOWER-1); 229 accumulator.msw = accumulator.midw = accumulator.lsw = 0;
260 230 /* Do the basic fixed point polynomial evaluation */
261 mul_Xsig_Xsig(&accumulator, &argSignif); 231 polynomial_Xsig(&accumulator, &Xsq, logterms, HIPOWER - 1);
262 shr_Xsig(&accumulator, 6 - adj); 232
263 233 mul_Xsig_Xsig(&accumulator, &argSignif);
264 mul32_Xsig(&arg_signif, leadterm); 234 shr_Xsig(&accumulator, 6 - adj);
265 add_two_Xsig(&accumulator, &arg_signif, &exponent); 235
266 236 mul32_Xsig(&arg_signif, leadterm);
267 *expon = exponent + 1; 237 add_two_Xsig(&accumulator, &arg_signif, &exponent);
268 accum_result->lsw = accumulator.lsw; 238
269 accum_result->midw = accumulator.midw; 239 *expon = exponent + 1;
270 accum_result->msw = accumulator.msw; 240 accum_result->lsw = accumulator.lsw;
241 accum_result->midw = accumulator.midw;
242 accum_result->msw = accumulator.msw;
271 243
272} 244}
generated by cgit v1.2.2 (git 2.25.0) at 2024-12-14 02:33:49 -0500