^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 1) // SPDX-License-Identifier: GPL-2.0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 2) /*---------------------------------------------------------------------------+
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 3) | poly_2xm1.c |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 4) | |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 5) | Function to compute 2^x-1 by a polynomial approximation. |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 6) | |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 7) | Copyright (C) 1992,1993,1994,1997 |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 8) | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 9) | E-mail billm@suburbia.net |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 10) | |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 11) | |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 12) +---------------------------------------------------------------------------*/
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 13)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 14) #include "exception.h"
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 15) #include "reg_constant.h"
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 16) #include "fpu_emu.h"
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 17) #include "fpu_system.h"
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 18) #include "control_w.h"
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 19) #include "poly.h"
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 20)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 21) #define HIPOWER 11
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 22) static const unsigned long long lterms[HIPOWER] = {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 23) 0x0000000000000000LL, /* This term done separately as 12 bytes */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 24) 0xf5fdeffc162c7543LL,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 25) 0x1c6b08d704a0bfa6LL,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 26) 0x0276556df749cc21LL,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 27) 0x002bb0ffcf14f6b8LL,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 28) 0x0002861225ef751cLL,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 29) 0x00001ffcbfcd5422LL,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 30) 0x00000162c005d5f1LL,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 31) 0x0000000da96ccb1bLL,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 32) 0x0000000078d1b897LL,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 33) 0x000000000422b029LL
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 34) };
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 35)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 36) static const Xsig hiterm = MK_XSIG(0xb17217f7, 0xd1cf79ab, 0xc8a39194);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 37)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 38) /* Four slices: 0.0 : 0.25 : 0.50 : 0.75 : 1.0,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 39) These numbers are 2^(1/4), 2^(1/2), and 2^(3/4)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 40) */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 41) static const Xsig shiftterm0 = MK_XSIG(0, 0, 0);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 42) static const Xsig shiftterm1 = MK_XSIG(0x9837f051, 0x8db8a96f, 0x46ad2318);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 43) static const Xsig shiftterm2 = MK_XSIG(0xb504f333, 0xf9de6484, 0x597d89b3);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 44) static const Xsig shiftterm3 = MK_XSIG(0xd744fcca, 0xd69d6af4, 0x39a68bb9);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 45)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 46) static const Xsig *shiftterm[] = { &shiftterm0, &shiftterm1,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 47) &shiftterm2, &shiftterm3
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 48) };
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 49)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 50) /*--- poly_2xm1() -----------------------------------------------------------+
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 51) | Requires st(0) which is TAG_Valid and < 1. |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 52) +---------------------------------------------------------------------------*/
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 53) int poly_2xm1(u_char sign, FPU_REG *arg, FPU_REG *result)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 54) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 55) long int exponent, shift;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 56) unsigned long long Xll;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 57) Xsig accumulator, Denom, argSignif;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 58) u_char tag;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 59)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 60) exponent = exponent16(arg);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 61)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 62) #ifdef PARANOID
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 63) if (exponent >= 0) { /* Don't want a |number| >= 1.0 */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 64) /* Number negative, too large, or not Valid. */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 65) EXCEPTION(EX_INTERNAL | 0x127);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 66) return 1;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 67) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 68) #endif /* PARANOID */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 69)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 70) argSignif.lsw = 0;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 71) XSIG_LL(argSignif) = Xll = significand(arg);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 72)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 73) if (exponent == -1) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 74) shift = (argSignif.msw & 0x40000000) ? 3 : 2;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 75) /* subtract 0.5 or 0.75 */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 76) exponent -= 2;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 77) XSIG_LL(argSignif) <<= 2;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 78) Xll <<= 2;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 79) } else if (exponent == -2) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 80) shift = 1;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 81) /* subtract 0.25 */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 82) exponent--;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 83) XSIG_LL(argSignif) <<= 1;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 84) Xll <<= 1;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 85) } else
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 86) shift = 0;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 87)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 88) if (exponent < -2) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 89) /* Shift the argument right by the required places. */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 90) if (FPU_shrx(&Xll, -2 - exponent) >= 0x80000000U)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 91) Xll++; /* round up */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 92) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 93)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 94) accumulator.lsw = accumulator.midw = accumulator.msw = 0;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 95) polynomial_Xsig(&accumulator, &Xll, lterms, HIPOWER - 1);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 96) mul_Xsig_Xsig(&accumulator, &argSignif);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 97) shr_Xsig(&accumulator, 3);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 98)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 99) mul_Xsig_Xsig(&argSignif, &hiterm); /* The leading term */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 100) add_two_Xsig(&accumulator, &argSignif, &exponent);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 101)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 102) if (shift) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 103) /* The argument is large, use the identity:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 104) f(x+a) = f(a) * (f(x) + 1) - 1;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 105) */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 106) shr_Xsig(&accumulator, -exponent);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 107) accumulator.msw |= 0x80000000; /* add 1.0 */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 108) mul_Xsig_Xsig(&accumulator, shiftterm[shift]);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 109) accumulator.msw &= 0x3fffffff; /* subtract 1.0 */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 110) exponent = 1;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 111) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 112)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 113) if (sign != SIGN_POS) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 114) /* The argument is negative, use the identity:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 115) f(-x) = -f(x) / (1 + f(x))
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 116) */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 117) Denom.lsw = accumulator.lsw;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 118) XSIG_LL(Denom) = XSIG_LL(accumulator);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 119) if (exponent < 0)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 120) shr_Xsig(&Denom, -exponent);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 121) else if (exponent > 0) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 122) /* exponent must be 1 here */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 123) XSIG_LL(Denom) <<= 1;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 124) if (Denom.lsw & 0x80000000)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 125) XSIG_LL(Denom) |= 1;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 126) (Denom.lsw) <<= 1;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 127) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 128) Denom.msw |= 0x80000000; /* add 1.0 */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 129) div_Xsig(&accumulator, &Denom, &accumulator);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 130) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 131)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 132) /* Convert to 64 bit signed-compatible */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 133) exponent += round_Xsig(&accumulator);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 134)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 135) result = &st(0);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 136) significand(result) = XSIG_LL(accumulator);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 137) setexponent16(result, exponent);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 138)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 139) tag = FPU_round(result, 1, 0, FULL_PRECISION, sign);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 140)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 141) setsign(result, sign);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 142) FPU_settag0(tag);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 143)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 144) return 0;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 145)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 146) }