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