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
+---------------------------------------------------------------------------*/
13
#include "exception.h"
14
#include "reg_constant.h"
16
#include "fpu_system.h"
17
#include "control_w.h"
20
static void log2_kernel(FPU_REG const *arg, u_char argsign,
21
Xsig * accum_result, long int *expon);
23
/*--- poly_l2() -------------------------------------------------------------+
24
| Base 2 logarithm by a polynomial approximation. |
25
+---------------------------------------------------------------------------*/
26
void poly_l2(FPU_REG *st0_ptr, FPU_REG *st1_ptr, u_char st1_sign)
28
long int exponent, expon, expon_expon;
29
Xsig accumulator, expon_accum, yaccum;
34
exponent = exponent16(st0_ptr);
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);
44
/* Treat as 1 <= st0_ptr < sqrt(2) */
45
x.sigh = st0_ptr->sigh - 0x80000000;
46
x.sigl = st0_ptr->sigl;
50
tag = FPU_normalize_nuo(&x);
52
if (tag == TAG_Zero) {
54
accumulator.msw = accumulator.midw = accumulator.lsw = 0;
56
log2_kernel(&x, argsign, &accumulator, &expon);
64
expon_accum.msw = exponent;
65
expon_accum.midw = expon_accum.lsw = 0;
67
expon_expon = 31 + norm_Xsig(&expon_accum);
68
shr_Xsig(&accumulator, expon_expon - expon);
71
negate_Xsig(&accumulator);
72
add_Xsig_Xsig(&accumulator, &expon_accum);
79
XSIG_LL(yaccum) = significand(st1_ptr);
80
mul_Xsig_Xsig(&accumulator, &yaccum);
82
expon_expon += round_Xsig(&accumulator);
84
if (accumulator.msw == 0) {
85
FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
89
significand(st1_ptr) = XSIG_LL(accumulator);
90
setexponent16(st1_ptr, expon_expon + exponent16(st1_ptr) + 1);
92
tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign ^ st1_sign);
95
set_precision_flag_up(); /* 80486 appears to always do this */
101
/*--- poly_l2p1() -----------------------------------------------------------+
102
| Base 2 logarithm by a polynomial approximation. |
104
+---------------------------------------------------------------------------*/
105
int poly_l2p1(u_char sign0, u_char sign1,
106
FPU_REG * st0_ptr, FPU_REG * st1_ptr, FPU_REG * dest)
110
Xsig accumulator, yaccum;
112
if (exponent16(st0_ptr) < 0) {
113
log2_kernel(st0_ptr, sign0, &accumulator, &exponent);
116
XSIG_LL(yaccum) = significand(st1_ptr);
117
mul_Xsig_Xsig(&accumulator, &yaccum);
119
exponent += round_Xsig(&accumulator);
121
exponent += exponent16(st1_ptr) + 1;
122
if (exponent < EXP_WAY_UNDER)
123
exponent = EXP_WAY_UNDER;
125
significand(dest) = XSIG_LL(accumulator);
126
setexponent16(dest, exponent);
128
tag = FPU_round(dest, 1, 0, FULL_PRECISION, sign0 ^ sign1);
131
if (tag == TAG_Valid)
132
set_precision_flag_up(); /* 80486 appears to always do this */
134
/* The magnitude of st0_ptr is far too large. */
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). */
141
if (arith_invalid(1) < 0)
143
#endif /* PECULIAR_486 */
146
/* 80486 appears to do this */
147
if (sign0 == SIGN_NEG)
148
set_precision_flag_down();
150
set_precision_flag_up();
153
if (exponent(dest) <= EXP_UNDER)
154
EXCEPTION(EX_Underflow);
162
static const unsigned long long logterms[HIPOWER] = {
163
0x2a8eca5705fc2ef0LL,
164
0xf6384ee1d01febceLL,
165
0x093bb62877cdf642LL,
166
0x006985d8a9ec439bLL,
167
0x0005212c4f55a9c8LL,
168
0x00004326a16927f0LL,
169
0x0000038d1d80a0e7LL,
170
0x0000003141cc80c6LL,
171
0x00000002b1668c9fLL,
175
static const unsigned long leadterm = 0xb8000000;
177
/*--- log2_kernel() ---------------------------------------------------------+
178
| Base 2 logarithm by a polynomial approximation. |
180
+---------------------------------------------------------------------------*/
181
static void log2_kernel(FPU_REG const *arg, u_char argsign, Xsig *accum_result,
184
long int exponent, adj;
185
unsigned long long Xsq;
186
Xsig accumulator, Numer, Denom, argSignif, arg_signif;
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);
196
shr_Xsig(&Denom, 1 - (1 + exponent));
198
if (Denom.msw & 0x80000000) {
199
div_Xsig(&Numer, &Denom, &argSignif);
202
/* Denom must be 1.0 */
203
argSignif.lsw = Numer.lsw;
204
argSignif.midw = Numer.midw;
205
argSignif.msw = Numer.msw;
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 */
216
#endif /* PECULIAR_486 */
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)
229
accumulator.msw = accumulator.midw = accumulator.lsw = 0;
230
/* Do the basic fixed point polynomial evaluation */
231
polynomial_Xsig(&accumulator, &Xsq, logterms, HIPOWER - 1);
233
mul_Xsig_Xsig(&accumulator, &argSignif);
234
shr_Xsig(&accumulator, 6 - adj);
236
mul32_Xsig(&arg_signif, leadterm);
237
add_two_Xsig(&accumulator, &arg_signif, &exponent);
239
*expon = exponent + 1;
240
accum_result->lsw = accumulator.lsw;
241
accum_result->midw = accumulator.midw;
242
accum_result->msw = accumulator.msw;