x86/math-emu/poly_l2.c

0001 // SPDX-License-Identifier: GPL-2.0
0002 /*---------------------------------------------------------------------------+
0003  |  poly_l2.c                                                                |
0004  |                                                                           |
0005  | Compute the base 2 log of a FPU_REG, using a polynomial approximation.    |
0006  |                                                                           |
0007  | Copyright (C) 1992,1993,1994,1997                                         |
0008  |                  W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
0009  |                  E-mail   billm@suburbia.net                              |
0010  |                                                                           |
0011  |                                                                           |
0012  +---------------------------------------------------------------------------*/
0013
0014 #include "exception.h"
0015 #include "reg_constant.h"
0016 #include "fpu_emu.h"
0017 #include "fpu_system.h"
0018 #include "control_w.h"
0019 #include "poly.h"
0020
0021 static void log2_kernel(FPU_REG const *arg, u_char argsign,
0022             Xsig * accum_result, long int *expon);
0023
0024 /*--- poly_l2() -------------------------------------------------------------+
0025  |   Base 2 logarithm by a polynomial approximation.                         |
0026  +---------------------------------------------------------------------------*/
0027 void poly_l2(FPU_REG *st0_ptr, FPU_REG *st1_ptr, u_char st1_sign)
0028 {
0029     long int exponent, expon, expon_expon;
0030     Xsig accumulator, expon_accum, yaccum;
0031     u_char sign, argsign;
0032     FPU_REG x;
0033     int tag;
0034
0035     exponent = exponent16(st0_ptr);
0036
0037     /* From st0_ptr, make a number > sqrt(2)/2 and < sqrt(2) */
0038     if (st0_ptr->sigh > (unsigned)0xb504f334) {
0039         /* Treat as  sqrt(2)/2 < st0_ptr < 1 */
0040         significand(&x) = -significand(st0_ptr);
0041         setexponent16(&x, -1);
0042         exponent++;
0043         argsign = SIGN_NEG;
0044     } else {
0045         /* Treat as  1 <= st0_ptr < sqrt(2) */
0046         x.sigh = st0_ptr->sigh - 0x80000000;
0047         x.sigl = st0_ptr->sigl;
0048         setexponent16(&x, 0);
0049         argsign = SIGN_POS;
0050     }
0051     tag = FPU_normalize_nuo(&x);
0052
0053     if (tag == TAG_Zero) {
0054         expon = 0;
0055         accumulator.msw = accumulator.midw = accumulator.lsw = 0;
0056     } else {
0057         log2_kernel(&x, argsign, &accumulator, &expon);
0058     }
0059
0060     if (exponent < 0) {
0061         sign = SIGN_NEG;
0062         exponent = -exponent;
0063     } else
0064         sign = SIGN_POS;
0065     expon_accum.msw = exponent;
0066     expon_accum.midw = expon_accum.lsw = 0;
0067     if (exponent) {
0068         expon_expon = 31 + norm_Xsig(&expon_accum);
0069         shr_Xsig(&accumulator, expon_expon - expon);
0070
0071         if (sign ^ argsign)
0072             negate_Xsig(&accumulator);
0073         add_Xsig_Xsig(&accumulator, &expon_accum);
0074     } else {
0075         expon_expon = expon;
0076         sign = argsign;
0077     }
0078
0079     yaccum.lsw = 0;
0080     XSIG_LL(yaccum) = significand(st1_ptr);
0081     mul_Xsig_Xsig(&accumulator, &yaccum);
0082
0083     expon_expon += round_Xsig(&accumulator);
0084
0085     if (accumulator.msw == 0) {
0086         FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
0087         return;
0088     }
0089
0090     significand(st1_ptr) = XSIG_LL(accumulator);
0091     setexponent16(st1_ptr, expon_expon + exponent16(st1_ptr) + 1);
0092
0093     tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign ^ st1_sign);
0094     FPU_settagi(1, tag);
0095
0096     set_precision_flag_up();    /* 80486 appears to always do this */
0097
0098     return;
0099
0100 }
0101
0102 /*--- poly_l2p1() -----------------------------------------------------------+
0103  |   Base 2 logarithm by a polynomial approximation.                         |
0104  |   log2(x+1)                                                               |
0105  +---------------------------------------------------------------------------*/
0106 int poly_l2p1(u_char sign0, u_char sign1,
0107           FPU_REG * st0_ptr, FPU_REG * st1_ptr, FPU_REG * dest)
0108 {
0109     u_char tag;
0110     long int exponent;
0111     Xsig accumulator, yaccum;
0112
0113     if (exponent16(st0_ptr) < 0) {
0114         log2_kernel(st0_ptr, sign0, &accumulator, &exponent);
0115
0116         yaccum.lsw = 0;
0117         XSIG_LL(yaccum) = significand(st1_ptr);
0118         mul_Xsig_Xsig(&accumulator, &yaccum);
0119
0120         exponent += round_Xsig(&accumulator);
0121
0122         exponent += exponent16(st1_ptr) + 1;
0123         if (exponent < EXP_WAY_UNDER)
0124             exponent = EXP_WAY_UNDER;
0125
0126         significand(dest) = XSIG_LL(accumulator);
0127         setexponent16(dest, exponent);
0128
0129         tag = FPU_round(dest, 1, 0, FULL_PRECISION, sign0 ^ sign1);
0130         FPU_settagi(1, tag);
0131
0132         if (tag == TAG_Valid)
0133             set_precision_flag_up();    /* 80486 appears to always do this */
0134     } else {
0135         /* The magnitude of st0_ptr is far too large. */
0136
0137         if (sign0 != SIGN_POS) {
0138             /* Trying to get the log of a negative number. */
0139 #ifdef PECULIAR_486     /* Stupid 80486 doesn't worry about log(negative). */
0140             changesign(st1_ptr);
0141 #else
0142             if (arith_invalid(1) < 0)
0143                 return 1;
0144 #endif /* PECULIAR_486 */
0145         }
0146
0147         /* 80486 appears to do this */
0148         if (sign0 == SIGN_NEG)
0149             set_precision_flag_down();
0150         else
0151             set_precision_flag_up();
0152     }
0153
0154     if (exponent(dest) <= EXP_UNDER)
0155         EXCEPTION(EX_Underflow);
0156
0157     return 0;
0158
0159 }
0160
0161 #undef HIPOWER
0162 #define HIPOWER 10
0163 static const unsigned long long logterms[HIPOWER] = {
0164     0x2a8eca5705fc2ef0LL,
0165     0xf6384ee1d01febceLL,
0166     0x093bb62877cdf642LL,
0167     0x006985d8a9ec439bLL,
0168     0x0005212c4f55a9c8LL,
0169     0x00004326a16927f0LL,
0170     0x0000038d1d80a0e7LL,
0171     0x0000003141cc80c6LL,
0172     0x00000002b1668c9fLL,
0173     0x000000002c7a46aaLL
0174 };
0175
0176 static const unsigned long leadterm = 0xb8000000;
0177
0178 /*--- log2_kernel() ---------------------------------------------------------+
0179  |   Base 2 logarithm by a polynomial approximation.                         |
0180  |   log2(x+1)                                                               |
0181  +---------------------------------------------------------------------------*/
0182 static void log2_kernel(FPU_REG const *arg, u_char argsign, Xsig *accum_result,
0183             long int *expon)
0184 {
0185     long int exponent, adj;
0186     unsigned long long Xsq;
0187     Xsig accumulator, Numer, Denom, argSignif, arg_signif;
0188
0189     exponent = exponent16(arg);
0190     Numer.lsw = Denom.lsw = 0;
0191     XSIG_LL(Numer) = XSIG_LL(Denom) = significand(arg);
0192     if (argsign == SIGN_POS) {
0193         shr_Xsig(&Denom, 2 - (1 + exponent));
0194         Denom.msw |= 0x80000000;
0195         div_Xsig(&Numer, &Denom, &argSignif);
0196     } else {
0197         shr_Xsig(&Denom, 1 - (1 + exponent));
0198         negate_Xsig(&Denom);
0199         if (Denom.msw & 0x80000000) {
0200             div_Xsig(&Numer, &Denom, &argSignif);
0201             exponent++;
0202         } else {
0203             /* Denom must be 1.0 */
0204             argSignif.lsw = Numer.lsw;
0205             argSignif.midw = Numer.midw;
0206             argSignif.msw = Numer.msw;
0207         }
0208     }
0209
0210 #ifndef PECULIAR_486
0211     /* Should check here that  |local_arg|  is within the valid range */
0212     if (exponent >= -2) {
0213         if ((exponent > -2) || (argSignif.msw > (unsigned)0xafb0ccc0)) {
0214             /* The argument is too large */
0215         }
0216     }
0217 #endif /* PECULIAR_486 */
0218
0219     arg_signif.lsw = argSignif.lsw;
0220     XSIG_LL(arg_signif) = XSIG_LL(argSignif);
0221     adj = norm_Xsig(&argSignif);
0222     accumulator.lsw = argSignif.lsw;
0223     XSIG_LL(accumulator) = XSIG_LL(argSignif);
0224     mul_Xsig_Xsig(&accumulator, &accumulator);
0225     shr_Xsig(&accumulator, 2 * (-1 - (1 + exponent + adj)));
0226     Xsq = XSIG_LL(accumulator);
0227     if (accumulator.lsw & 0x80000000)
0228         Xsq++;
0229
0230     accumulator.msw = accumulator.midw = accumulator.lsw = 0;
0231     /* Do the basic fixed point polynomial evaluation */
0232     polynomial_Xsig(&accumulator, &Xsq, logterms, HIPOWER - 1);
0233
0234     mul_Xsig_Xsig(&accumulator, &argSignif);
0235     shr_Xsig(&accumulator, 6 - adj);
0236
0237     mul32_Xsig(&arg_signif, leadterm);
0238     add_two_Xsig(&accumulator, &arg_signif, &exponent);
0239
0240     *expon = exponent + 1;
0241     accum_result->lsw = accumulator.lsw;
0242     accum_result->midw = accumulator.midw;
0243     accum_result->msw = accumulator.msw;
0244
0245 }