1 /*---------------------------------------------------------------------------+
4 | Compute the base 2 log of a FPU_REG, using a polynomial approximation. |
6 | Copyright (C) 1992,1993,1994,1997 |
7 | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
8 | E-mail billm@suburbia.net |
11 +---------------------------------------------------------------------------*/
14 #include "exception.h"
15 #include "reg_constant.h"
17 #include "fpu_system.h"
18 #include "control_w.h"
22 static void log2_kernel(FPU_REG
const *arg
, u_char argsign
,
23 Xsig
*accum_result
, long int *expon
);
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
)
31 long int exponent
, expon
, expon_expon
;
32 Xsig accumulator
, expon_accum
, yaccum
;
37 exponent
= exponent16(st0_ptr
);
39 /* From st0_ptr, make a number > sqrt(2)/2 and < sqrt(2) */
40 if ( st0_ptr
->sigh
> (unsigned)0xb504f334 )
42 /* Treat as sqrt(2)/2 < st0_ptr < 1 */
43 significand(&x
) = - significand(st0_ptr
);
44 setexponent16(&x
, -1);
50 /* Treat as 1 <= st0_ptr < sqrt(2) */
51 x
.sigh
= st0_ptr
->sigh
- 0x80000000;
52 x
.sigl
= st0_ptr
->sigl
;
56 tag
= FPU_normalize_nuo(&x
);
58 if ( tag
== TAG_Zero
)
61 accumulator
.msw
= accumulator
.midw
= accumulator
.lsw
= 0;
65 log2_kernel(&x
, argsign
, &accumulator
, &expon
);
75 expon_accum
.msw
= exponent
; expon_accum
.midw
= expon_accum
.lsw
= 0;
78 expon_expon
= 31 + norm_Xsig(&expon_accum
);
79 shr_Xsig(&accumulator
, expon_expon
- expon
);
82 negate_Xsig(&accumulator
);
83 add_Xsig_Xsig(&accumulator
, &expon_accum
);
91 yaccum
.lsw
= 0; XSIG_LL(yaccum
) = significand(st1_ptr
);
92 mul_Xsig_Xsig(&accumulator
, &yaccum
);
94 expon_expon
+= round_Xsig(&accumulator
);
96 if ( accumulator
.msw
== 0 )
98 FPU_copy_to_reg1(&CONST_Z
, TAG_Zero
);
102 significand(st1_ptr
) = XSIG_LL(accumulator
);
103 setexponent16(st1_ptr
, expon_expon
+ exponent16(st1_ptr
) + 1);
105 tag
= FPU_round(st1_ptr
, 1, 0, FULL_PRECISION
, sign
^ st1_sign
);
108 set_precision_flag_up(); /* 80486 appears to always do this */
115 /*--- poly_l2p1() -----------------------------------------------------------+
116 | Base 2 logarithm by a polynomial approximation. |
118 +---------------------------------------------------------------------------*/
119 int poly_l2p1(u_char sign0
, u_char sign1
,
120 FPU_REG
*st0_ptr
, FPU_REG
*st1_ptr
, FPU_REG
*dest
)
124 Xsig accumulator
, yaccum
;
126 if ( exponent16(st0_ptr
) < 0 )
128 log2_kernel(st0_ptr
, sign0
, &accumulator
, &exponent
);
131 XSIG_LL(yaccum
) = significand(st1_ptr
);
132 mul_Xsig_Xsig(&accumulator
, &yaccum
);
134 exponent
+= round_Xsig(&accumulator
);
136 exponent
+= exponent16(st1_ptr
) + 1;
137 if ( exponent
< EXP_WAY_UNDER
) exponent
= EXP_WAY_UNDER
;
139 significand(dest
) = XSIG_LL(accumulator
);
140 setexponent16(dest
, exponent
);
142 tag
= FPU_round(dest
, 1, 0, FULL_PRECISION
, sign0
^ sign1
);
145 if ( tag
== TAG_Valid
)
146 set_precision_flag_up(); /* 80486 appears to always do this */
150 /* The magnitude of st0_ptr is far too large. */
152 if ( sign0
!= SIGN_POS
)
154 /* Trying to get the log of a negative number. */
155 #ifdef PECULIAR_486 /* Stupid 80486 doesn't worry about log(negative). */
158 if ( arith_invalid(1) < 0 )
163 /* 80486 appears to do this */
164 if ( sign0
== SIGN_NEG
)
165 set_precision_flag_down();
167 set_precision_flag_up();
170 if ( exponent(dest
) <= EXP_UNDER
)
171 EXCEPTION(EX_Underflow
);
182 static const unsigned long long logterms
[HIPOWER
] =
184 0x2a8eca5705fc2ef0LL
,
185 0xf6384ee1d01febceLL
,
186 0x093bb62877cdf642LL
,
187 0x006985d8a9ec439bLL
,
188 0x0005212c4f55a9c8LL
,
189 0x00004326a16927f0LL
,
190 0x0000038d1d80a0e7LL
,
191 0x0000003141cc80c6LL
,
192 0x00000002b1668c9fLL
,
196 static const unsigned long leadterm
= 0xb8000000;
199 /*--- log2_kernel() ---------------------------------------------------------+
200 | Base 2 logarithm by a polynomial approximation. |
202 +---------------------------------------------------------------------------*/
203 static void log2_kernel(FPU_REG
const *arg
, u_char argsign
, Xsig
*accum_result
,
206 long int exponent
, adj
;
207 unsigned long long Xsq
;
208 Xsig accumulator
, Numer
, Denom
, argSignif
, arg_signif
;
210 exponent
= exponent16(arg
);
211 Numer
.lsw
= Denom
.lsw
= 0;
212 XSIG_LL(Numer
) = XSIG_LL(Denom
) = significand(arg
);
213 if ( argsign
== SIGN_POS
)
215 shr_Xsig(&Denom
, 2 - (1 + exponent
));
216 Denom
.msw
|= 0x80000000;
217 div_Xsig(&Numer
, &Denom
, &argSignif
);
221 shr_Xsig(&Denom
, 1 - (1 + exponent
));
223 if ( Denom
.msw
& 0x80000000 )
225 div_Xsig(&Numer
, &Denom
, &argSignif
);
230 /* Denom must be 1.0 */
231 argSignif
.lsw
= Numer
.lsw
; argSignif
.midw
= Numer
.midw
;
232 argSignif
.msw
= Numer
.msw
;
237 /* Should check here that |local_arg| is within the valid range */
238 if ( exponent
>= -2 )
240 if ( (exponent
> -2) ||
241 (argSignif
.msw
> (unsigned)0xafb0ccc0) )
243 /* The argument is too large */
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 )
257 accumulator
.msw
= accumulator
.midw
= accumulator
.lsw
= 0;
258 /* Do the basic fixed point polynomial evaluation */
259 polynomial_Xsig(&accumulator
, &Xsq
, logterms
, HIPOWER
-1);
261 mul_Xsig_Xsig(&accumulator
, &argSignif
);
262 shr_Xsig(&accumulator
, 6 - adj
);
264 mul32_Xsig(&arg_signif
, leadterm
);
265 add_two_Xsig(&accumulator
, &arg_signif
, &exponent
);
267 *expon
= exponent
+ 1;
268 accum_result
->lsw
= accumulator
.lsw
;
269 accum_result
->midw
= accumulator
.midw
;
270 accum_result
->msw
= accumulator
.msw
;