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 20 static 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 +---------------------------------------------------------------------------*/ 26 void 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 +---------------------------------------------------------------------------*/ 105 int 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 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, 172 0x000000002c7a46aaLL 173 }; 174 175 static const unsigned long leadterm = 0xb8000000; 176 177 /*--- log2_kernel() ---------------------------------------------------------+ 178 | Base 2 logarithm by a polynomial approximation. | 179 | log2(x+1) | 180 +---------------------------------------------------------------------------*/ 181 static 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 } 245