1da957e11SThomas Gleixner /*---------------------------------------------------------------------------+ 2da957e11SThomas Gleixner | poly_atan.c | 3da957e11SThomas Gleixner | | 4da957e11SThomas Gleixner | Compute the arctan of a FPU_REG, using a polynomial approximation. | 5da957e11SThomas Gleixner | | 6da957e11SThomas Gleixner | Copyright (C) 1992,1993,1994,1997 | 7da957e11SThomas Gleixner | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia | 8da957e11SThomas Gleixner | E-mail billm@suburbia.net | 9da957e11SThomas Gleixner | | 10da957e11SThomas Gleixner | | 11da957e11SThomas Gleixner +---------------------------------------------------------------------------*/ 12da957e11SThomas Gleixner 13da957e11SThomas Gleixner #include "exception.h" 14da957e11SThomas Gleixner #include "reg_constant.h" 15da957e11SThomas Gleixner #include "fpu_emu.h" 16da957e11SThomas Gleixner #include "fpu_system.h" 17da957e11SThomas Gleixner #include "status_w.h" 18da957e11SThomas Gleixner #include "control_w.h" 19da957e11SThomas Gleixner #include "poly.h" 20da957e11SThomas Gleixner 21da957e11SThomas Gleixner 22da957e11SThomas Gleixner #define HIPOWERon 6 /* odd poly, negative terms */ 23da957e11SThomas Gleixner static const unsigned long long oddnegterms[HIPOWERon] = 24da957e11SThomas Gleixner { 25da957e11SThomas Gleixner 0x0000000000000000LL, /* Dummy (not for - 1.0) */ 26da957e11SThomas Gleixner 0x015328437f756467LL, 27da957e11SThomas Gleixner 0x0005dda27b73dec6LL, 28da957e11SThomas Gleixner 0x0000226bf2bfb91aLL, 29da957e11SThomas Gleixner 0x000000ccc439c5f7LL, 30da957e11SThomas Gleixner 0x0000000355438407LL 31da957e11SThomas Gleixner } ; 32da957e11SThomas Gleixner 33da957e11SThomas Gleixner #define HIPOWERop 6 /* odd poly, positive terms */ 34da957e11SThomas Gleixner static const unsigned long long oddplterms[HIPOWERop] = 35da957e11SThomas Gleixner { 36da957e11SThomas Gleixner /* 0xaaaaaaaaaaaaaaabLL, transferred to fixedpterm[] */ 37da957e11SThomas Gleixner 0x0db55a71875c9ac2LL, 38da957e11SThomas Gleixner 0x0029fce2d67880b0LL, 39da957e11SThomas Gleixner 0x0000dfd3908b4596LL, 40da957e11SThomas Gleixner 0x00000550fd61dab4LL, 41da957e11SThomas Gleixner 0x0000001c9422b3f9LL, 42da957e11SThomas Gleixner 0x000000003e3301e1LL 43da957e11SThomas Gleixner }; 44da957e11SThomas Gleixner 45da957e11SThomas Gleixner static const unsigned long long denomterm = 0xebd9b842c5c53a0eLL; 46da957e11SThomas Gleixner 47da957e11SThomas Gleixner static const Xsig fixedpterm = MK_XSIG(0xaaaaaaaa, 0xaaaaaaaa, 0xaaaaaaaa); 48da957e11SThomas Gleixner 49da957e11SThomas Gleixner static const Xsig pi_signif = MK_XSIG(0xc90fdaa2, 0x2168c234, 0xc4c6628b); 50da957e11SThomas Gleixner 51da957e11SThomas Gleixner 52da957e11SThomas Gleixner /*--- poly_atan() -----------------------------------------------------------+ 53da957e11SThomas Gleixner | | 54da957e11SThomas Gleixner +---------------------------------------------------------------------------*/ 55da957e11SThomas Gleixner void poly_atan(FPU_REG *st0_ptr, u_char st0_tag, 56da957e11SThomas Gleixner FPU_REG *st1_ptr, u_char st1_tag) 57da957e11SThomas Gleixner { 58da957e11SThomas Gleixner u_char transformed, inverted, 59da957e11SThomas Gleixner sign1, sign2; 60da957e11SThomas Gleixner int exponent; 61da957e11SThomas Gleixner long int dummy_exp; 62da957e11SThomas Gleixner Xsig accumulator, Numer, Denom, accumulatore, argSignif, 63da957e11SThomas Gleixner argSq, argSqSq; 64da957e11SThomas Gleixner u_char tag; 65da957e11SThomas Gleixner 66da957e11SThomas Gleixner sign1 = getsign(st0_ptr); 67da957e11SThomas Gleixner sign2 = getsign(st1_ptr); 68da957e11SThomas Gleixner if ( st0_tag == TAG_Valid ) 69da957e11SThomas Gleixner { 70da957e11SThomas Gleixner exponent = exponent(st0_ptr); 71da957e11SThomas Gleixner } 72da957e11SThomas Gleixner else 73da957e11SThomas Gleixner { 74da957e11SThomas Gleixner /* This gives non-compatible stack contents... */ 75da957e11SThomas Gleixner FPU_to_exp16(st0_ptr, st0_ptr); 76da957e11SThomas Gleixner exponent = exponent16(st0_ptr); 77da957e11SThomas Gleixner } 78da957e11SThomas Gleixner if ( st1_tag == TAG_Valid ) 79da957e11SThomas Gleixner { 80da957e11SThomas Gleixner exponent -= exponent(st1_ptr); 81da957e11SThomas Gleixner } 82da957e11SThomas Gleixner else 83da957e11SThomas Gleixner { 84da957e11SThomas Gleixner /* This gives non-compatible stack contents... */ 85da957e11SThomas Gleixner FPU_to_exp16(st1_ptr, st1_ptr); 86da957e11SThomas Gleixner exponent -= exponent16(st1_ptr); 87da957e11SThomas Gleixner } 88da957e11SThomas Gleixner 89da957e11SThomas Gleixner if ( (exponent < 0) || ((exponent == 0) && 90da957e11SThomas Gleixner ((st0_ptr->sigh < st1_ptr->sigh) || 91da957e11SThomas Gleixner ((st0_ptr->sigh == st1_ptr->sigh) && 92da957e11SThomas Gleixner (st0_ptr->sigl < st1_ptr->sigl))) ) ) 93da957e11SThomas Gleixner { 94da957e11SThomas Gleixner inverted = 1; 95da957e11SThomas Gleixner Numer.lsw = Denom.lsw = 0; 96da957e11SThomas Gleixner XSIG_LL(Numer) = significand(st0_ptr); 97da957e11SThomas Gleixner XSIG_LL(Denom) = significand(st1_ptr); 98da957e11SThomas Gleixner } 99da957e11SThomas Gleixner else 100da957e11SThomas Gleixner { 101da957e11SThomas Gleixner inverted = 0; 102da957e11SThomas Gleixner exponent = -exponent; 103da957e11SThomas Gleixner Numer.lsw = Denom.lsw = 0; 104da957e11SThomas Gleixner XSIG_LL(Numer) = significand(st1_ptr); 105da957e11SThomas Gleixner XSIG_LL(Denom) = significand(st0_ptr); 106da957e11SThomas Gleixner } 107da957e11SThomas Gleixner div_Xsig(&Numer, &Denom, &argSignif); 108da957e11SThomas Gleixner exponent += norm_Xsig(&argSignif); 109da957e11SThomas Gleixner 110da957e11SThomas Gleixner if ( (exponent >= -1) 111da957e11SThomas Gleixner || ((exponent == -2) && (argSignif.msw > 0xd413ccd0)) ) 112da957e11SThomas Gleixner { 113da957e11SThomas Gleixner /* The argument is greater than sqrt(2)-1 (=0.414213562...) */ 114da957e11SThomas Gleixner /* Convert the argument by an identity for atan */ 115da957e11SThomas Gleixner transformed = 1; 116da957e11SThomas Gleixner 117da957e11SThomas Gleixner if ( exponent >= 0 ) 118da957e11SThomas Gleixner { 119da957e11SThomas Gleixner #ifdef PARANOID 120da957e11SThomas Gleixner if ( !( (exponent == 0) && 121da957e11SThomas Gleixner (argSignif.lsw == 0) && (argSignif.midw == 0) && 122da957e11SThomas Gleixner (argSignif.msw == 0x80000000) ) ) 123da957e11SThomas Gleixner { 124da957e11SThomas Gleixner EXCEPTION(EX_INTERNAL|0x104); /* There must be a logic error */ 125da957e11SThomas Gleixner return; 126da957e11SThomas Gleixner } 127da957e11SThomas Gleixner #endif /* PARANOID */ 128da957e11SThomas Gleixner argSignif.msw = 0; /* Make the transformed arg -> 0.0 */ 129da957e11SThomas Gleixner } 130da957e11SThomas Gleixner else 131da957e11SThomas Gleixner { 132da957e11SThomas Gleixner Numer.lsw = Denom.lsw = argSignif.lsw; 133da957e11SThomas Gleixner XSIG_LL(Numer) = XSIG_LL(Denom) = XSIG_LL(argSignif); 134da957e11SThomas Gleixner 135da957e11SThomas Gleixner if ( exponent < -1 ) 136da957e11SThomas Gleixner shr_Xsig(&Numer, -1-exponent); 137da957e11SThomas Gleixner negate_Xsig(&Numer); 138da957e11SThomas Gleixner 139da957e11SThomas Gleixner shr_Xsig(&Denom, -exponent); 140da957e11SThomas Gleixner Denom.msw |= 0x80000000; 141da957e11SThomas Gleixner 142da957e11SThomas Gleixner div_Xsig(&Numer, &Denom, &argSignif); 143da957e11SThomas Gleixner 144da957e11SThomas Gleixner exponent = -1 + norm_Xsig(&argSignif); 145da957e11SThomas Gleixner } 146da957e11SThomas Gleixner } 147da957e11SThomas Gleixner else 148da957e11SThomas Gleixner { 149da957e11SThomas Gleixner transformed = 0; 150da957e11SThomas Gleixner } 151da957e11SThomas Gleixner 152da957e11SThomas Gleixner argSq.lsw = argSignif.lsw; argSq.midw = argSignif.midw; 153da957e11SThomas Gleixner argSq.msw = argSignif.msw; 154da957e11SThomas Gleixner mul_Xsig_Xsig(&argSq, &argSq); 155da957e11SThomas Gleixner 156da957e11SThomas Gleixner argSqSq.lsw = argSq.lsw; argSqSq.midw = argSq.midw; argSqSq.msw = argSq.msw; 157da957e11SThomas Gleixner mul_Xsig_Xsig(&argSqSq, &argSqSq); 158da957e11SThomas Gleixner 159da957e11SThomas Gleixner accumulatore.lsw = argSq.lsw; 160da957e11SThomas Gleixner XSIG_LL(accumulatore) = XSIG_LL(argSq); 161da957e11SThomas Gleixner 162da957e11SThomas Gleixner shr_Xsig(&argSq, 2*(-1-exponent-1)); 163da957e11SThomas Gleixner shr_Xsig(&argSqSq, 4*(-1-exponent-1)); 164da957e11SThomas Gleixner 165da957e11SThomas Gleixner /* Now have argSq etc with binary point at the left 166da957e11SThomas Gleixner .1xxxxxxxx */ 167da957e11SThomas Gleixner 168da957e11SThomas Gleixner /* Do the basic fixed point polynomial evaluation */ 169da957e11SThomas Gleixner accumulator.msw = accumulator.midw = accumulator.lsw = 0; 170da957e11SThomas Gleixner polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq), 171da957e11SThomas Gleixner oddplterms, HIPOWERop-1); 172da957e11SThomas Gleixner mul64_Xsig(&accumulator, &XSIG_LL(argSq)); 173da957e11SThomas Gleixner negate_Xsig(&accumulator); 174da957e11SThomas Gleixner polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq), oddnegterms, HIPOWERon-1); 175da957e11SThomas Gleixner negate_Xsig(&accumulator); 176da957e11SThomas Gleixner add_two_Xsig(&accumulator, &fixedpterm, &dummy_exp); 177da957e11SThomas Gleixner 178da957e11SThomas Gleixner mul64_Xsig(&accumulatore, &denomterm); 179da957e11SThomas Gleixner shr_Xsig(&accumulatore, 1 + 2*(-1-exponent)); 180da957e11SThomas Gleixner accumulatore.msw |= 0x80000000; 181da957e11SThomas Gleixner 182da957e11SThomas Gleixner div_Xsig(&accumulator, &accumulatore, &accumulator); 183da957e11SThomas Gleixner 184da957e11SThomas Gleixner mul_Xsig_Xsig(&accumulator, &argSignif); 185da957e11SThomas Gleixner mul_Xsig_Xsig(&accumulator, &argSq); 186da957e11SThomas Gleixner 187da957e11SThomas Gleixner shr_Xsig(&accumulator, 3); 188da957e11SThomas Gleixner negate_Xsig(&accumulator); 189da957e11SThomas Gleixner add_Xsig_Xsig(&accumulator, &argSignif); 190da957e11SThomas Gleixner 191da957e11SThomas Gleixner if ( transformed ) 192da957e11SThomas Gleixner { 193da957e11SThomas Gleixner /* compute pi/4 - accumulator */ 194da957e11SThomas Gleixner shr_Xsig(&accumulator, -1-exponent); 195da957e11SThomas Gleixner negate_Xsig(&accumulator); 196da957e11SThomas Gleixner add_Xsig_Xsig(&accumulator, &pi_signif); 197da957e11SThomas Gleixner exponent = -1; 198da957e11SThomas Gleixner } 199da957e11SThomas Gleixner 200da957e11SThomas Gleixner if ( inverted ) 201da957e11SThomas Gleixner { 202da957e11SThomas Gleixner /* compute pi/2 - accumulator */ 203da957e11SThomas Gleixner shr_Xsig(&accumulator, -exponent); 204da957e11SThomas Gleixner negate_Xsig(&accumulator); 205da957e11SThomas Gleixner add_Xsig_Xsig(&accumulator, &pi_signif); 206da957e11SThomas Gleixner exponent = 0; 207da957e11SThomas Gleixner } 208da957e11SThomas Gleixner 209da957e11SThomas Gleixner if ( sign1 ) 210da957e11SThomas Gleixner { 211da957e11SThomas Gleixner /* compute pi - accumulator */ 212da957e11SThomas Gleixner shr_Xsig(&accumulator, 1 - exponent); 213da957e11SThomas Gleixner negate_Xsig(&accumulator); 214da957e11SThomas Gleixner add_Xsig_Xsig(&accumulator, &pi_signif); 215da957e11SThomas Gleixner exponent = 1; 216da957e11SThomas Gleixner } 217da957e11SThomas Gleixner 218da957e11SThomas Gleixner exponent += round_Xsig(&accumulator); 219da957e11SThomas Gleixner 220da957e11SThomas Gleixner significand(st1_ptr) = XSIG_LL(accumulator); 221da957e11SThomas Gleixner setexponent16(st1_ptr, exponent); 222da957e11SThomas Gleixner 223da957e11SThomas Gleixner tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign2); 224da957e11SThomas Gleixner FPU_settagi(1, tag); 225da957e11SThomas Gleixner 226da957e11SThomas Gleixner set_precision_flag_up(); /* We do not really know if up or down, 227da957e11SThomas Gleixner use this as the default. */ 228da957e11SThomas Gleixner 229da957e11SThomas Gleixner } 230