11541Srgrimes// SPDX-License-Identifier: GPL-2.0 21541Srgrimes/*---------------------------------------------------------------------------+ 31541Srgrimes | poly_2xm1.c | 41541Srgrimes | | 51541Srgrimes | Function to compute 2^x-1 by a polynomial approximation. | 61541Srgrimes | | 71541Srgrimes | Copyright (C) 1992,1993,1994,1997 | 81541Srgrimes | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia | 91541Srgrimes | E-mail billm@suburbia.net | 101541Srgrimes | | 111541Srgrimes | | 121541Srgrimes +---------------------------------------------------------------------------*/ 131541Srgrimes 141541Srgrimes#include "exception.h" 151541Srgrimes#include "reg_constant.h" 161541Srgrimes#include "fpu_emu.h" 171541Srgrimes#include "fpu_system.h" 181541Srgrimes#include "control_w.h" 191541Srgrimes#include "poly.h" 201541Srgrimes 211541Srgrimes#define HIPOWER 11 221541Srgrimesstatic const unsigned long long lterms[HIPOWER] = { 231541Srgrimes 0x0000000000000000LL, /* This term done separately as 12 bytes */ 241541Srgrimes 0xf5fdeffc162c7543LL, 251541Srgrimes 0x1c6b08d704a0bfa6LL, 261541Srgrimes 0x0276556df749cc21LL, 271541Srgrimes 0x002bb0ffcf14f6b8LL, 281541Srgrimes 0x0002861225ef751cLL, 291541Srgrimes 0x00001ffcbfcd5422LL, 301541Srgrimes 0x00000162c005d5f1LL, 311541Srgrimes 0x0000000da96ccb1bLL, 321541Srgrimes 0x0000000078d1b897LL, 3322521Sdyson 0x000000000422b029LL 3450477Speter}; 351541Srgrimes 361541Srgrimesstatic const Xsig hiterm = MK_XSIG(0xb17217f7, 0xd1cf79ab, 0xc8a39194); 372165Spaul 382165Spaul/* Four slices: 0.0 : 0.25 : 0.50 : 0.75 : 1.0, 392165Spaul These numbers are 2^(1/4), 2^(1/2), and 2^(3/4) 401541Srgrimes */ 411541Srgrimesstatic const Xsig shiftterm0 = MK_XSIG(0, 0, 0); 4269564Salfredstatic const Xsig shiftterm1 = MK_XSIG(0x9837f051, 0x8db8a96f, 0x46ad2318); 43106582Smuxstatic const Xsig shiftterm2 = MK_XSIG(0xb504f333, 0xf9de6484, 0x597d89b3); 44101848Srwatsonstatic const Xsig shiftterm3 = MK_XSIG(0xd744fcca, 0xd69d6af4, 0x39a68bb9); 4576166Smarkm 4676166Smarkmstatic const Xsig *shiftterm[] = { &shiftterm0, &shiftterm1, 4769564Salfred &shiftterm2, &shiftterm3 481541Srgrimes}; 4996755Strhodes 501541Srgrimes/*--- poly_2xm1() -----------------------------------------------------------+ 511541Srgrimes | Requires st(0) which is TAG_Valid and < 1. | 521541Srgrimes +---------------------------------------------------------------------------*/ 531541Srgrimesint poly_2xm1(u_char sign, FPU_REG *arg, FPU_REG *result) 541541Srgrimes{ 551541Srgrimes long int exponent, shift; 561541Srgrimes unsigned long long Xll; 571541Srgrimes Xsig accumulator, Denom, argSignif; 581541Srgrimes u_char tag; 591541Srgrimes 601541Srgrimes exponent = exponent16(arg); 611541Srgrimes 621541Srgrimes#ifdef PARANOID 631541Srgrimes if (exponent >= 0) { /* Don't want a |number| >= 1.0 */ 6496755Strhodes /* Number negative, too large, or not Valid. */ 651541Srgrimes EXCEPTION(EX_INTERNAL | 0x127); 661541Srgrimes return 1; 67106576Smux } 68106576Smux#endif /* PARANOID */ 691541Srgrimes 70106582Smux argSignif.lsw = 0; 71106582Smux XSIG_LL(argSignif) = Xll = significand(arg); 721541Srgrimes 7318006Sdg if (exponent == -1) { 7496755Strhodes shift = (argSignif.msw & 0x40000000) ? 3 : 2; 751541Srgrimes /* subtract 0.5 or 0.75 */ 7696755Strhodes exponent -= 2; 771541Srgrimes XSIG_LL(argSignif) <<= 2; 781541Srgrimes Xll <<= 2; 7996755Strhodes } else if (exponent == -2) { 801541Srgrimes shift = 1; 8196755Strhodes /* subtract 0.25 */ 8218006Sdg exponent--; 8341169Sbde XSIG_LL(argSignif) <<= 1; 8431132Sjulian Xll <<= 1; 85106582Smux } else 86106582Smux shift = 0; 8722521Sdyson 881541Srgrimes if (exponent < -2) { 89106582Smux /* Shift the argument right by the required places. */ 90106582Smux if (FPU_shrx(&Xll, -2 - exponent) >= 0x80000000U) 9153975Smckusick Xll++; /* round up */ 921541Srgrimes } 9353975Smckusick 94106582Smux accumulator.lsw = accumulator.midw = accumulator.msw = 0; 95106582Smux polynomial_Xsig(&accumulator, &Xll, lterms, HIPOWER - 1); 96106582Smux mul_Xsig_Xsig(&accumulator, &argSignif); 97106582Smux shr_Xsig(&accumulator, 3); 98106582Smux 99106582Smux mul_Xsig_Xsig(&argSignif, &hiterm); /* The leading term */ 1001541Srgrimes add_two_Xsig(&accumulator, &argSignif, &exponent); 1011541Srgrimes 10269564Salfred if (shift) { 103106582Smux /* The argument is large, use the identity: 104106582Smux f(x+a) = f(a) * (f(x) + 1) - 1; 105106582Smux */ 106106582Smux shr_Xsig(&accumulator, -exponent); 107106582Smux accumulator.msw |= 0x80000000; /* add 1.0 */ 108106582Smux mul_Xsig_Xsig(&accumulator, shiftterm[shift]); 109106582Smux accumulator.msw &= 0x3fffffff; /* subtract 1.0 */ 110106582Smux exponent = 1; 111106582Smux } 112106582Smux 113106582Smux if (sign != SIGN_POS) { 1141541Srgrimes /* The argument is negative, use the identity: 11596755Strhodes f(-x) = -f(x) / (1 + f(x)) 11696755Strhodes */ 1171541Srgrimes Denom.lsw = accumulator.lsw; 11886037Sdillon XSIG_LL(Denom) = XSIG_LL(accumulator); 11986037Sdillon if (exponent < 0) 12086037Sdillon shr_Xsig(&Denom, -exponent); 12186037Sdillon else if (exponent > 0) { 12286037Sdillon /* exponent must be 1 here */ 12386037Sdillon XSIG_LL(Denom) <<= 1; 12486037Sdillon if (Denom.lsw & 0x80000000) 1251541Srgrimes XSIG_LL(Denom) |= 1; 1261541Srgrimes (Denom.lsw) <<= 1; 12760938Sjake } 1281541Srgrimes Denom.msw |= 0x80000000; /* add 1.0 */ 12910027Sdg div_Xsig(&accumulator, &Denom, &accumulator); 1301541Srgrimes } 13134266Sjulian 13285339Sdillon /* Convert to 64 bit signed-compatible */ 13386037Sdillon exponent += round_Xsig(&accumulator); 13422521Sdyson 13562553Smckusick result = &st(0); 13631144Sjulian significand(result) = XSIG_LL(accumulator); 13797188Smux setexponent16(result, exponent); 13897188Smux 13931132Sjulian tag = FPU_round(result, 1, 0, FULL_PRECISION, sign); 1401541Srgrimes 1411541Srgrimes setsign(result, sign); 142102088Sphk FPU_settag0(tag); 1431541Srgrimes 14410358Sjulian return 0; 145106582Smux 14675934Sphk} 147100985Srwatson