1 // SPDX-License-Identifier: GPL-2.0 2 /*---------------------------------------------------------------------------+ 3 | poly_l2.c | 4 | | 5 | Compute the base 2 log of a FPU_REG, using a polynomial approximation. | 6 | | 7 | Copyright (C) 1992,1993,1994,1997 | 8 | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia | 9 | E-mail billm@suburbia.net | 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 static void log2_kernel(FPU_REG const *arg, u_char argsign, 22 Xsig * accum_result, long int *expon); 23 24 /*--- poly_l2() -------------------------------------------------------------+ 25 | Base 2 logarithm by a polynomial approximation. | 26 +---------------------------------------------------------------------------*/ 27 void poly_l2(FPU_REG *st0_ptr, FPU_REG *st1_ptr, u_char st1_sign) 28 { 29 long int exponent, expon, expon_expon; 30 Xsig accumulator, expon_accum, yaccum; 31 u_char sign, argsign; 32 FPU_REG x; 33 int tag; 34 35 exponent = exponent16(st0_ptr); 36 37 /* From st0_ptr, make a number > sqrt(2)/2 and < sqrt(2) */ 38 if (st0_ptr->sigh > (unsigned)0xb504f334) { 39 /* Treat as sqrt(2)/2 < st0_ptr < 1 */ 40 significand(&x) = -significand(st0_ptr); 41 setexponent16(&x, -1); 42 exponent++; 43 argsign = SIGN_NEG; 44 } else { 45 /* Treat as 1 <= st0_ptr < sqrt(2) */ 46 x.sigh = st0_ptr->sigh - 0x80000000; 47 x.sigl = st0_ptr->sigl; 48 setexponent16(&x, 0); 49 argsign = SIGN_POS; 50 } 51 tag = FPU_normalize_nuo(&x); 52 53 if (tag == TAG_Zero) { 54 expon = 0; 55 accumulator.msw = accumulator.midw = accumulator.lsw = 0; 56 } else { 57 log2_kernel(&x, argsign, &accumulator, &expon); 58 } 59 60 if (exponent < 0) { 61 sign = SIGN_NEG; 62 exponent = -exponent; 63 } else 64 sign = SIGN_POS; 65 expon_accum.msw = exponent; 66 expon_accum.midw = expon_accum.lsw = 0; 67 if (exponent) { 68 expon_expon = 31 + norm_Xsig(&expon_accum); 69 shr_Xsig(&accumulator, expon_expon - expon); 70 71 if (sign ^ argsign) 72 negate_Xsig(&accumulator); 73 add_Xsig_Xsig(&accumulator, &expon_accum); 74 } else { 75 expon_expon = expon; 76 sign = argsign; 77 } 78 79 yaccum.lsw = 0; 80 XSIG_LL(yaccum) = significand(st1_ptr); 81 mul_Xsig_Xsig(&accumulator, &yaccum); 82 83 expon_expon += round_Xsig(&accumulator); 84 85 if (accumulator.msw == 0) { 86 FPU_copy_to_reg1(&CONST_Z, TAG_Zero); 87 return; 88 } 89 90 significand(st1_ptr) = XSIG_LL(accumulator); 91 setexponent16(st1_ptr, expon_expon + exponent16(st1_ptr) + 1); 92 93 tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign ^ st1_sign); 94 FPU_settagi(1, tag); 95 96 set_precision_flag_up(); /* 80486 appears to always do this */ 97 98 return; 99 100 } 101 102 /*--- poly_l2p1() -----------------------------------------------------------+ 103 | Base 2 logarithm by a polynomial approximation. | 104 | log2(x+1) | 105 +---------------------------------------------------------------------------*/ 106 int poly_l2p1(u_char sign0, u_char sign1, 107 FPU_REG * st0_ptr, FPU_REG * st1_ptr, FPU_REG * dest) 108 { 109 u_char tag; 110 long int exponent; 111 Xsig accumulator, yaccum; 112 113 if (exponent16(st0_ptr) < 0) { 114 log2_kernel(st0_ptr, sign0, &accumulator, &exponent); 115 116 yaccum.lsw = 0; 117 XSIG_LL(yaccum) = significand(st1_ptr); 118 mul_Xsig_Xsig(&accumulator, &yaccum); 119 120 exponent += round_Xsig(&accumulator); 121 122 exponent += exponent16(st1_ptr) + 1; 123 if (exponent < EXP_WAY_UNDER) 124 exponent = EXP_WAY_UNDER; 125 126 significand(dest) = XSIG_LL(accumulator); 127 setexponent16(dest, exponent); 128 129 tag = FPU_round(dest, 1, 0, FULL_PRECISION, sign0 ^ sign1); 130 FPU_settagi(1, tag); 131 132 if (tag == TAG_Valid) 133 set_precision_flag_up(); /* 80486 appears to always do this */ 134 } else { 135 /* The magnitude of st0_ptr is far too large. */ 136 137 if (sign0 != SIGN_POS) { 138 /* Trying to get the log of a negative number. */ 139 #ifdef PECULIAR_486 /* Stupid 80486 doesn't worry about log(negative). */ 140 changesign(st1_ptr); 141 #else 142 if (arith_invalid(1) < 0) 143 return 1; 144 #endif /* PECULIAR_486 */ 145 } 146 147 /* 80486 appears to do this */ 148 if (sign0 == SIGN_NEG) 149 set_precision_flag_down(); 150 else 151 set_precision_flag_up(); 152 } 153 154 if (exponent(dest) <= EXP_UNDER) 155 EXCEPTION(EX_Underflow); 156 157 return 0; 158 159 } 160 161 #undef HIPOWER 162 #define HIPOWER 10 163 static const unsigned long long logterms[HIPOWER] = { 164 0x2a8eca5705fc2ef0LL, 165 0xf6384ee1d01febceLL, 166 0x093bb62877cdf642LL, 167 0x006985d8a9ec439bLL, 168 0x0005212c4f55a9c8LL, 169 0x00004326a16927f0LL, 170 0x0000038d1d80a0e7LL, 171 0x0000003141cc80c6LL, 172 0x00000002b1668c9fLL, 173 0x000000002c7a46aaLL 174 }; 175 176 static const unsigned long leadterm = 0xb8000000; 177 178 /*--- log2_kernel() ---------------------------------------------------------+ 179 | Base 2 logarithm by a polynomial approximation. | 180 | log2(x+1) | 181 +---------------------------------------------------------------------------*/ 182 static void log2_kernel(FPU_REG const *arg, u_char argsign, Xsig *accum_result, 183 long int *expon) 184 { 185 long int exponent, adj; 186 unsigned long long Xsq; 187 Xsig accumulator, Numer, Denom, argSignif, arg_signif; 188 189 exponent = exponent16(arg); 190 Numer.lsw = Denom.lsw = 0; 191 XSIG_LL(Numer) = XSIG_LL(Denom) = significand(arg); 192 if (argsign == SIGN_POS) { 193 shr_Xsig(&Denom, 2 - (1 + exponent)); 194 Denom.msw |= 0x80000000; 195 div_Xsig(&Numer, &Denom, &argSignif); 196 } else { 197 shr_Xsig(&Denom, 1 - (1 + exponent)); 198 negate_Xsig(&Denom); 199 if (Denom.msw & 0x80000000) { 200 div_Xsig(&Numer, &Denom, &argSignif); 201 exponent++; 202 } else { 203 /* Denom must be 1.0 */ 204 argSignif.lsw = Numer.lsw; 205 argSignif.midw = Numer.midw; 206 argSignif.msw = Numer.msw; 207 } 208 } 209 210 #ifndef PECULIAR_486 211 /* Should check here that |local_arg| is within the valid range */ 212 if (exponent >= -2) { 213 if ((exponent > -2) || (argSignif.msw > (unsigned)0xafb0ccc0)) { 214 /* The argument is too large */ 215 } 216 } 217 #endif /* PECULIAR_486 */ 218 219 arg_signif.lsw = argSignif.lsw; 220 XSIG_LL(arg_signif) = XSIG_LL(argSignif); 221 adj = norm_Xsig(&argSignif); 222 accumulator.lsw = argSignif.lsw; 223 XSIG_LL(accumulator) = XSIG_LL(argSignif); 224 mul_Xsig_Xsig(&accumulator, &accumulator); 225 shr_Xsig(&accumulator, 2 * (-1 - (1 + exponent + adj))); 226 Xsq = XSIG_LL(accumulator); 227 if (accumulator.lsw & 0x80000000) 228 Xsq++; 229 230 accumulator.msw = accumulator.midw = accumulator.lsw = 0; 231 /* Do the basic fixed point polynomial evaluation */ 232 polynomial_Xsig(&accumulator, &Xsq, logterms, HIPOWER - 1); 233 234 mul_Xsig_Xsig(&accumulator, &argSignif); 235 shr_Xsig(&accumulator, 6 - adj); 236 237 mul32_Xsig(&arg_signif, leadterm); 238 add_two_Xsig(&accumulator, &arg_signif, &exponent); 239 240 *expon = exponent + 1; 241 accum_result->lsw = accumulator.lsw; 242 accum_result->midw = accumulator.midw; 243 accum_result->msw = accumulator.msw; 244 245 } 246