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