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