From ace9429bb58fd418f0c81d4c2835699bddf6bde6 Mon Sep 17 00:00:00 2001 From: Daniel Baumann Date: Thu, 11 Apr 2024 10:27:49 +0200 Subject: Adding upstream version 6.6.15. Signed-off-by: Daniel Baumann --- arch/x86/math-emu/poly_l2.c | 245 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 245 insertions(+) create mode 100644 arch/x86/math-emu/poly_l2.c (limited to 'arch/x86/math-emu/poly_l2.c') diff --git a/arch/x86/math-emu/poly_l2.c b/arch/x86/math-emu/poly_l2.c new file mode 100644 index 0000000000..98b6949bb8 --- /dev/null +++ b/arch/x86/math-emu/poly_l2.c @@ -0,0 +1,245 @@ +// SPDX-License-Identifier: GPL-2.0 +/*---------------------------------------------------------------------------+ + | poly_l2.c | + | | + | Compute the base 2 log of a FPU_REG, using a polynomial approximation. | + | | + | Copyright (C) 1992,1993,1994,1997 | + | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia | + | E-mail billm@suburbia.net | + | | + | | + +---------------------------------------------------------------------------*/ + +#include "exception.h" +#include "reg_constant.h" +#include "fpu_emu.h" +#include "fpu_system.h" +#include "control_w.h" +#include "poly.h" + +static void log2_kernel(FPU_REG const *arg, u_char argsign, + Xsig * accum_result, long int *expon); + +/*--- poly_l2() -------------------------------------------------------------+ + | Base 2 logarithm by a polynomial approximation. | + +---------------------------------------------------------------------------*/ +void poly_l2(FPU_REG *st0_ptr, FPU_REG *st1_ptr, u_char st1_sign) +{ + long int exponent, expon, expon_expon; + Xsig accumulator, expon_accum, yaccum; + u_char sign, argsign; + FPU_REG x; + int tag; + + exponent = exponent16(st0_ptr); + + /* From st0_ptr, make a number > sqrt(2)/2 and < sqrt(2) */ + if (st0_ptr->sigh > (unsigned)0xb504f334) { + /* Treat as sqrt(2)/2 < st0_ptr < 1 */ + significand(&x) = -significand(st0_ptr); + setexponent16(&x, -1); + exponent++; + argsign = SIGN_NEG; + } else { + /* Treat as 1 <= st0_ptr < sqrt(2) */ + x.sigh = st0_ptr->sigh - 0x80000000; + x.sigl = st0_ptr->sigl; + setexponent16(&x, 0); + argsign = SIGN_POS; + } + tag = FPU_normalize_nuo(&x); + + if (tag == TAG_Zero) { + expon = 0; + accumulator.msw = accumulator.midw = accumulator.lsw = 0; + } else { + log2_kernel(&x, argsign, &accumulator, &expon); + } + + if (exponent < 0) { + sign = SIGN_NEG; + exponent = -exponent; + } else + sign = SIGN_POS; + expon_accum.msw = exponent; + expon_accum.midw = expon_accum.lsw = 0; + if (exponent) { + expon_expon = 31 + norm_Xsig(&expon_accum); + shr_Xsig(&accumulator, expon_expon - expon); + + if (sign ^ argsign) + negate_Xsig(&accumulator); + add_Xsig_Xsig(&accumulator, &expon_accum); + } else { + expon_expon = expon; + sign = argsign; + } + + yaccum.lsw = 0; + XSIG_LL(yaccum) = significand(st1_ptr); + mul_Xsig_Xsig(&accumulator, &yaccum); + + expon_expon += round_Xsig(&accumulator); + + if (accumulator.msw == 0) { + FPU_copy_to_reg1(&CONST_Z, TAG_Zero); + return; + } + + significand(st1_ptr) = XSIG_LL(accumulator); + setexponent16(st1_ptr, expon_expon + exponent16(st1_ptr) + 1); + + tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign ^ st1_sign); + FPU_settagi(1, tag); + + set_precision_flag_up(); /* 80486 appears to always do this */ + + return; + +} + +/*--- poly_l2p1() -----------------------------------------------------------+ + | Base 2 logarithm by a polynomial approximation. | + | log2(x+1) | + +---------------------------------------------------------------------------*/ +int poly_l2p1(u_char sign0, u_char sign1, + FPU_REG * st0_ptr, FPU_REG * st1_ptr, FPU_REG * dest) +{ + u_char tag; + long int exponent; + Xsig accumulator, yaccum; + + if (exponent16(st0_ptr) < 0) { + log2_kernel(st0_ptr, sign0, &accumulator, &exponent); + + yaccum.lsw = 0; + XSIG_LL(yaccum) = significand(st1_ptr); + mul_Xsig_Xsig(&accumulator, &yaccum); + + exponent += round_Xsig(&accumulator); + + exponent += exponent16(st1_ptr) + 1; + if (exponent < EXP_WAY_UNDER) + exponent = EXP_WAY_UNDER; + + significand(dest) = XSIG_LL(accumulator); + setexponent16(dest, exponent); + + tag = FPU_round(dest, 1, 0, FULL_PRECISION, sign0 ^ sign1); + FPU_settagi(1, tag); + + if (tag == TAG_Valid) + set_precision_flag_up(); /* 80486 appears to always do this */ + } else { + /* The magnitude of st0_ptr is far too large. */ + + if (sign0 != SIGN_POS) { + /* Trying to get the log of a negative number. */ +#ifdef PECULIAR_486 /* Stupid 80486 doesn't worry about log(negative). */ + changesign(st1_ptr); +#else + if (arith_invalid(1) < 0) + return 1; +#endif /* PECULIAR_486 */ + } + + /* 80486 appears to do this */ + if (sign0 == SIGN_NEG) + set_precision_flag_down(); + else + set_precision_flag_up(); + } + + if (exponent(dest) <= EXP_UNDER) + EXCEPTION(EX_Underflow); + + return 0; + +} + +#undef HIPOWER +#define HIPOWER 10 +static const unsigned long long logterms[HIPOWER] = { + 0x2a8eca5705fc2ef0LL, + 0xf6384ee1d01febceLL, + 0x093bb62877cdf642LL, + 0x006985d8a9ec439bLL, + 0x0005212c4f55a9c8LL, + 0x00004326a16927f0LL, + 0x0000038d1d80a0e7LL, + 0x0000003141cc80c6LL, + 0x00000002b1668c9fLL, + 0x000000002c7a46aaLL +}; + +static const unsigned long leadterm = 0xb8000000; + +/*--- log2_kernel() ---------------------------------------------------------+ + | Base 2 logarithm by a polynomial approximation. | + | log2(x+1) | + +---------------------------------------------------------------------------*/ +static void log2_kernel(FPU_REG const *arg, u_char argsign, Xsig *accum_result, + long int *expon) +{ + long int exponent, adj; + unsigned long long Xsq; + Xsig accumulator, Numer, Denom, argSignif, arg_signif; + + exponent = exponent16(arg); + Numer.lsw = Denom.lsw = 0; + XSIG_LL(Numer) = XSIG_LL(Denom) = significand(arg); + if (argsign == SIGN_POS) { + shr_Xsig(&Denom, 2 - (1 + exponent)); + Denom.msw |= 0x80000000; + div_Xsig(&Numer, &Denom, &argSignif); + } else { + shr_Xsig(&Denom, 1 - (1 + exponent)); + negate_Xsig(&Denom); + if (Denom.msw & 0x80000000) { + div_Xsig(&Numer, &Denom, &argSignif); + exponent++; + } else { + /* Denom must be 1.0 */ + argSignif.lsw = Numer.lsw; + argSignif.midw = Numer.midw; + argSignif.msw = Numer.msw; + } + } + +#ifndef PECULIAR_486 + /* Should check here that |local_arg| is within the valid range */ + if (exponent >= -2) { + if ((exponent > -2) || (argSignif.msw > (unsigned)0xafb0ccc0)) { + /* The argument is too large */ + } + } +#endif /* PECULIAR_486 */ + + arg_signif.lsw = argSignif.lsw; + XSIG_LL(arg_signif) = XSIG_LL(argSignif); + adj = norm_Xsig(&argSignif); + accumulator.lsw = argSignif.lsw; + XSIG_LL(accumulator) = XSIG_LL(argSignif); + mul_Xsig_Xsig(&accumulator, &accumulator); + shr_Xsig(&accumulator, 2 * (-1 - (1 + exponent + adj))); + Xsq = XSIG_LL(accumulator); + if (accumulator.lsw & 0x80000000) + Xsq++; + + accumulator.msw = accumulator.midw = accumulator.lsw = 0; + /* Do the basic fixed point polynomial evaluation */ + polynomial_Xsig(&accumulator, &Xsq, logterms, HIPOWER - 1); + + mul_Xsig_Xsig(&accumulator, &argSignif); + shr_Xsig(&accumulator, 6 - adj); + + mul32_Xsig(&arg_signif, leadterm); + add_two_Xsig(&accumulator, &arg_signif, &exponent); + + *expon = exponent + 1; + accum_result->lsw = accumulator.lsw; + accum_result->midw = accumulator.midw; + accum_result->msw = accumulator.msw; + +} -- cgit v1.2.3