1181880Sjhb// SPDX-License-Identifier: GPL-2.0 2181880Sjhb/*---------------------------------------------------------------------------+ 3181880Sjhb | poly_atan.c | 4181880Sjhb | | 5181880Sjhb | Compute the arctan of a FPU_REG, using a polynomial approximation. | 6181880Sjhb | | 7181880Sjhb | Copyright (C) 1992,1993,1994,1997 | 8181880Sjhb | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia | 9181880Sjhb | E-mail billm@suburbia.net | 10181880Sjhb | | 11181880Sjhb | | 12181880Sjhb +---------------------------------------------------------------------------*/ 13181880Sjhb 14181880Sjhb#include "exception.h" 15181880Sjhb#include "reg_constant.h" 16181880Sjhb#include "fpu_emu.h" 17181880Sjhb#include "fpu_system.h" 18181880Sjhb#include "status_w.h" 19181880Sjhb#include "control_w.h" 20181880Sjhb#include "poly.h" 21181880Sjhb 22181880Sjhb#define HIPOWERon 6 /* odd poly, negative terms */ 23181880Sjhbstatic const unsigned long long oddnegterms[HIPOWERon] = { 24181880Sjhb 0x0000000000000000LL, /* Dummy (not for - 1.0) */ 25181880Sjhb 0x015328437f756467LL, 26181880Sjhb 0x0005dda27b73dec6LL, 27181880Sjhb 0x0000226bf2bfb91aLL, 28181880Sjhb 0x000000ccc439c5f7LL, 29181880Sjhb 0x0000000355438407LL 30181880Sjhb}; 31181880Sjhb 32181880Sjhb#define HIPOWERop 6 /* odd poly, positive terms */ 33181880Sjhbstatic const unsigned long long oddplterms[HIPOWERop] = { 34181880Sjhb/* 0xaaaaaaaaaaaaaaabLL, transferred to fixedpterm[] */ 35181880Sjhb 0x0db55a71875c9ac2LL, 36181880Sjhb 0x0029fce2d67880b0LL, 37181880Sjhb 0x0000dfd3908b4596LL, 38181880Sjhb 0x00000550fd61dab4LL, 39181880Sjhb 0x0000001c9422b3f9LL, 40181880Sjhb 0x000000003e3301e1LL 41181880Sjhb}; 42194186Sed 43181880Sjhbstatic const unsigned long long denomterm = 0xebd9b842c5c53a0eLL; 44181880Sjhb 45181880Sjhbstatic const Xsig fixedpterm = MK_XSIG(0xaaaaaaaa, 0xaaaaaaaa, 0xaaaaaaaa); 46181880Sjhb 47181880Sjhbstatic const Xsig pi_signif = MK_XSIG(0xc90fdaa2, 0x2168c234, 0xc4c6628b); 48181880Sjhb 49181880Sjhb/*--- poly_atan() -----------------------------------------------------------+ 50181880Sjhb | | 51181880Sjhb +---------------------------------------------------------------------------*/ 52181880Sjhbvoid poly_atan(FPU_REG *st0_ptr, u_char st0_tag, 53181880Sjhb FPU_REG *st1_ptr, u_char st1_tag) 54181880Sjhb{ 55181880Sjhb u_char transformed, inverted, sign1, sign2; 56181880Sjhb int exponent; 57181880Sjhb long int dummy_exp; 58181880Sjhb Xsig accumulator, Numer, Denom, accumulatore, argSignif, argSq, argSqSq; 59181880Sjhb u_char tag; 60181880Sjhb 61181880Sjhb sign1 = getsign(st0_ptr); 62181880Sjhb sign2 = getsign(st1_ptr); 63181880Sjhb if (st0_tag == TAG_Valid) { 64181880Sjhb exponent = exponent(st0_ptr); 65181880Sjhb } else { 66181880Sjhb /* This gives non-compatible stack contents... */ 67181880Sjhb FPU_to_exp16(st0_ptr, st0_ptr); 68181880Sjhb exponent = exponent16(st0_ptr); 69181880Sjhb } 70181880Sjhb if (st1_tag == TAG_Valid) { 71181880Sjhb exponent -= exponent(st1_ptr); 72181880Sjhb } else { 73181880Sjhb /* This gives non-compatible stack contents... */ 74181880Sjhb FPU_to_exp16(st1_ptr, st1_ptr); 75181880Sjhb exponent -= exponent16(st1_ptr); 76181880Sjhb } 77181880Sjhb 78181880Sjhb if ((exponent < 0) || ((exponent == 0) && 79181880Sjhb ((st0_ptr->sigh < st1_ptr->sigh) || 80181880Sjhb ((st0_ptr->sigh == st1_ptr->sigh) && 81181880Sjhb (st0_ptr->sigl < st1_ptr->sigl))))) { 82181880Sjhb inverted = 1; 83181880Sjhb Numer.lsw = Denom.lsw = 0; 84181880Sjhb XSIG_LL(Numer) = significand(st0_ptr); 85181880Sjhb XSIG_LL(Denom) = significand(st1_ptr); 86181880Sjhb } else { 87181880Sjhb inverted = 0; 88181880Sjhb exponent = -exponent; 89181880Sjhb Numer.lsw = Denom.lsw = 0; 90181880Sjhb XSIG_LL(Numer) = significand(st1_ptr); 91181880Sjhb XSIG_LL(Denom) = significand(st0_ptr); 92181880Sjhb } 93181880Sjhb div_Xsig(&Numer, &Denom, &argSignif); 94181880Sjhb exponent += norm_Xsig(&argSignif); 95181880Sjhb 96181880Sjhb if ((exponent >= -1) 97181880Sjhb || ((exponent == -2) && (argSignif.msw > 0xd413ccd0))) { 98181880Sjhb /* The argument is greater than sqrt(2)-1 (=0.414213562...) */ 99181880Sjhb /* Convert the argument by an identity for atan */ 100181880Sjhb transformed = 1; 101181880Sjhb 102181880Sjhb if (exponent >= 0) { 103181880Sjhb#ifdef PARANOID 104181880Sjhb if (!((exponent == 0) && 105181880Sjhb (argSignif.lsw == 0) && (argSignif.midw == 0) && 106181880Sjhb (argSignif.msw == 0x80000000))) { 107181880Sjhb EXCEPTION(EX_INTERNAL | 0x104); /* There must be a logic error */ 108181880Sjhb return; 109181880Sjhb } 110181880Sjhb#endif /* PARANOID */ 111181880Sjhb argSignif.msw = 0; /* Make the transformed arg -> 0.0 */ 112181880Sjhb } else { 113181880Sjhb Numer.lsw = Denom.lsw = argSignif.lsw; 114181880Sjhb XSIG_LL(Numer) = XSIG_LL(Denom) = XSIG_LL(argSignif); 115181880Sjhb 116181880Sjhb if (exponent < -1) 117181880Sjhb shr_Xsig(&Numer, -1 - exponent); 118181880Sjhb negate_Xsig(&Numer); 119181880Sjhb 120181880Sjhb shr_Xsig(&Denom, -exponent); 121181880Sjhb Denom.msw |= 0x80000000; 122181880Sjhb 123181880Sjhb div_Xsig(&Numer, &Denom, &argSignif); 124181880Sjhb 125181880Sjhb exponent = -1 + norm_Xsig(&argSignif); 126181880Sjhb } 127181880Sjhb } else { 128181880Sjhb transformed = 0; 129181880Sjhb } 130181880Sjhb 131181880Sjhb argSq.lsw = argSignif.lsw; 132181880Sjhb argSq.midw = argSignif.midw; 133181880Sjhb argSq.msw = argSignif.msw; 134 mul_Xsig_Xsig(&argSq, &argSq); 135 136 argSqSq.lsw = argSq.lsw; 137 argSqSq.midw = argSq.midw; 138 argSqSq.msw = argSq.msw; 139 mul_Xsig_Xsig(&argSqSq, &argSqSq); 140 141 accumulatore.lsw = argSq.lsw; 142 XSIG_LL(accumulatore) = XSIG_LL(argSq); 143 144 shr_Xsig(&argSq, 2 * (-1 - exponent - 1)); 145 shr_Xsig(&argSqSq, 4 * (-1 - exponent - 1)); 146 147 /* Now have argSq etc with binary point at the left 148 .1xxxxxxxx */ 149 150 /* Do the basic fixed point polynomial evaluation */ 151 accumulator.msw = accumulator.midw = accumulator.lsw = 0; 152 polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq), 153 oddplterms, HIPOWERop - 1); 154 mul64_Xsig(&accumulator, &XSIG_LL(argSq)); 155 negate_Xsig(&accumulator); 156 polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq), oddnegterms, 157 HIPOWERon - 1); 158 negate_Xsig(&accumulator); 159 add_two_Xsig(&accumulator, &fixedpterm, &dummy_exp); 160 161 mul64_Xsig(&accumulatore, &denomterm); 162 shr_Xsig(&accumulatore, 1 + 2 * (-1 - exponent)); 163 accumulatore.msw |= 0x80000000; 164 165 div_Xsig(&accumulator, &accumulatore, &accumulator); 166 167 mul_Xsig_Xsig(&accumulator, &argSignif); 168 mul_Xsig_Xsig(&accumulator, &argSq); 169 170 shr_Xsig(&accumulator, 3); 171 negate_Xsig(&accumulator); 172 add_Xsig_Xsig(&accumulator, &argSignif); 173 174 if (transformed) { 175 /* compute pi/4 - accumulator */ 176 shr_Xsig(&accumulator, -1 - exponent); 177 negate_Xsig(&accumulator); 178 add_Xsig_Xsig(&accumulator, &pi_signif); 179 exponent = -1; 180 } 181 182 if (inverted) { 183 /* compute pi/2 - accumulator */ 184 shr_Xsig(&accumulator, -exponent); 185 negate_Xsig(&accumulator); 186 add_Xsig_Xsig(&accumulator, &pi_signif); 187 exponent = 0; 188 } 189 190 if (sign1) { 191 /* compute pi - accumulator */ 192 shr_Xsig(&accumulator, 1 - exponent); 193 negate_Xsig(&accumulator); 194 add_Xsig_Xsig(&accumulator, &pi_signif); 195 exponent = 1; 196 } 197 198 exponent += round_Xsig(&accumulator); 199 200 significand(st1_ptr) = XSIG_LL(accumulator); 201 setexponent16(st1_ptr, exponent); 202 203 tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign2); 204 FPU_settagi(1, tag); 205 206 set_precision_flag_up(); /* We do not really know if up or down, 207 use this as the default. */ 208 209} 210