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