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