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