1 /*---------------------------------------------------------------------------+ 2 | poly_tan.c | 3 | | 4 | Compute the tan of a FPU_REG, using a polynomial approximation. | 5 | | 6 | Copyright (C) 1992,1993,1994,1997,1999 | 7 | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, | 8 | Australia. E-mail billm@melbpc.org.au | 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 "control_w.h" 18 #include "poly.h" 19 20 21 #define HiPOWERop 3 /* odd poly, positive terms */ 22 static const unsigned long long oddplterm[HiPOWERop] = 23 { 24 0x0000000000000000LL, 25 0x0051a1cf08fca228LL, 26 0x0000000071284ff7LL 27 }; 28 29 #define HiPOWERon 2 /* odd poly, negative terms */ 30 static const unsigned long long oddnegterm[HiPOWERon] = 31 { 32 0x1291a9a184244e80LL, 33 0x0000583245819c21LL 34 }; 35 36 #define HiPOWERep 2 /* even poly, positive terms */ 37 static const unsigned long long evenplterm[HiPOWERep] = 38 { 39 0x0e848884b539e888LL, 40 0x00003c7f18b887daLL 41 }; 42 43 #define HiPOWERen 2 /* even poly, negative terms */ 44 static const unsigned long long evennegterm[HiPOWERen] = 45 { 46 0xf1f0200fd51569ccLL, 47 0x003afb46105c4432LL 48 }; 49 50 static const unsigned long long twothirds = 0xaaaaaaaaaaaaaaabLL; 51 52 53 /*--- poly_tan() ------------------------------------------------------------+ 54 | | 55 +---------------------------------------------------------------------------*/ 56 void poly_tan(FPU_REG *st0_ptr) 57 { 58 long int exponent; 59 int invert; 60 Xsig argSq, argSqSq, accumulatoro, accumulatore, accum, 61 argSignif, fix_up; 62 unsigned long adj; 63 64 exponent = exponent(st0_ptr); 65 66 #ifdef PARANOID 67 if ( signnegative(st0_ptr) ) /* Can't hack a number < 0.0 */ 68 { arith_invalid(0); return; } /* Need a positive number */ 69 #endif /* PARANOID */ 70 71 /* Split the problem into two domains, smaller and larger than pi/4 */ 72 if ( (exponent == 0) || ((exponent == -1) && (st0_ptr->sigh > 0xc90fdaa2)) ) 73 { 74 /* The argument is greater than (approx) pi/4 */ 75 invert = 1; 76 accum.lsw = 0; 77 XSIG_LL(accum) = significand(st0_ptr); 78 79 if ( exponent == 0 ) 80 { 81 /* The argument is >= 1.0 */ 82 /* Put the binary point at the left. */ 83 XSIG_LL(accum) <<= 1; 84 } 85 /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */ 86 XSIG_LL(accum) = 0x921fb54442d18469LL - XSIG_LL(accum); 87 /* This is a special case which arises due to rounding. */ 88 if ( XSIG_LL(accum) == 0xffffffffffffffffLL ) 89 { 90 FPU_settag0(TAG_Valid); 91 significand(st0_ptr) = 0x8a51e04daabda360LL; 92 setexponent16(st0_ptr, (0x41 + EXTENDED_Ebias) | SIGN_Negative); 93 return; 94 } 95 96 argSignif.lsw = accum.lsw; 97 XSIG_LL(argSignif) = XSIG_LL(accum); 98 exponent = -1 + norm_Xsig(&argSignif); 99 } 100 else 101 { 102 invert = 0; 103 argSignif.lsw = 0; 104 XSIG_LL(accum) = XSIG_LL(argSignif) = significand(st0_ptr); 105 106 if ( exponent < -1 ) 107 { 108 /* shift the argument right by the required places */ 109 if ( FPU_shrx(&XSIG_LL(accum), -1-exponent) >= 0x80000000U ) 110 XSIG_LL(accum) ++; /* round up */ 111 } 112 } 113 114 XSIG_LL(argSq) = XSIG_LL(accum); argSq.lsw = accum.lsw; 115 mul_Xsig_Xsig(&argSq, &argSq); 116 XSIG_LL(argSqSq) = XSIG_LL(argSq); argSqSq.lsw = argSq.lsw; 117 mul_Xsig_Xsig(&argSqSq, &argSqSq); 118 119 /* Compute the negative terms for the numerator polynomial */ 120 accumulatoro.msw = accumulatoro.midw = accumulatoro.lsw = 0; 121 polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddnegterm, HiPOWERon-1); 122 mul_Xsig_Xsig(&accumulatoro, &argSq); 123 negate_Xsig(&accumulatoro); 124 /* Add the positive terms */ 125 polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddplterm, HiPOWERop-1); 126 127 128 /* Compute the positive terms for the denominator polynomial */ 129 accumulatore.msw = accumulatore.midw = accumulatore.lsw = 0; 130 polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evenplterm, HiPOWERep-1); 131 mul_Xsig_Xsig(&accumulatore, &argSq); 132 negate_Xsig(&accumulatore); 133 /* Add the negative terms */ 134 polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evennegterm, HiPOWERen-1); 135 /* Multiply by arg^2 */ 136 mul64_Xsig(&accumulatore, &XSIG_LL(argSignif)); 137 mul64_Xsig(&accumulatore, &XSIG_LL(argSignif)); 138 /* de-normalize and divide by 2 */ 139 shr_Xsig(&accumulatore, -2*(1+exponent) + 1); 140 negate_Xsig(&accumulatore); /* This does 1 - accumulator */ 141 142 /* Now find the ratio. */ 143 if ( accumulatore.msw == 0 ) 144 { 145 /* accumulatoro must contain 1.0 here, (actually, 0) but it 146 really doesn't matter what value we use because it will 147 have negligible effect in later calculations 148 */ 149 XSIG_LL(accum) = 0x8000000000000000LL; 150 accum.lsw = 0; 151 } 152 else 153 { 154 div_Xsig(&accumulatoro, &accumulatore, &accum); 155 } 156 157 /* Multiply by 1/3 * arg^3 */ 158 mul64_Xsig(&accum, &XSIG_LL(argSignif)); 159 mul64_Xsig(&accum, &XSIG_LL(argSignif)); 160 mul64_Xsig(&accum, &XSIG_LL(argSignif)); 161 mul64_Xsig(&accum, &twothirds); 162 shr_Xsig(&accum, -2*(exponent+1)); 163 164 /* tan(arg) = arg + accum */ 165 add_two_Xsig(&accum, &argSignif, &exponent); 166 167 if ( invert ) 168 { 169 /* We now have the value of tan(pi_2 - arg) where pi_2 is an 170 approximation for pi/2 171 */ 172 /* The next step is to fix the answer to compensate for the 173 error due to the approximation used for pi/2 174 */ 175 176 /* This is (approx) delta, the error in our approx for pi/2 177 (see above). It has an exponent of -65 178 */ 179 XSIG_LL(fix_up) = 0x898cc51701b839a2LL; 180 fix_up.lsw = 0; 181 182 if ( exponent == 0 ) 183 adj = 0xffffffff; /* We want approx 1.0 here, but 184 this is close enough. */ 185 else if ( exponent > -30 ) 186 { 187 adj = accum.msw >> -(exponent+1); /* tan */ 188 adj = mul_32_32(adj, adj); /* tan^2 */ 189 } 190 else 191 adj = 0; 192 adj = mul_32_32(0x898cc517, adj); /* delta * tan^2 */ 193 194 fix_up.msw += adj; 195 if ( !(fix_up.msw & 0x80000000) ) /* did fix_up overflow ? */ 196 { 197 /* Yes, we need to add an msb */ 198 shr_Xsig(&fix_up, 1); 199 fix_up.msw |= 0x80000000; 200 shr_Xsig(&fix_up, 64 + exponent); 201 } 202 else 203 shr_Xsig(&fix_up, 65 + exponent); 204 205 add_two_Xsig(&accum, &fix_up, &exponent); 206 207 /* accum now contains tan(pi/2 - arg). 208 Use tan(arg) = 1.0 / tan(pi/2 - arg) 209 */ 210 accumulatoro.lsw = accumulatoro.midw = 0; 211 accumulatoro.msw = 0x80000000; 212 div_Xsig(&accumulatoro, &accum, &accum); 213 exponent = - exponent - 1; 214 } 215 216 /* Transfer the result */ 217 round_Xsig(&accum); 218 FPU_settag0(TAG_Valid); 219 significand(st0_ptr) = XSIG_LL(accum); 220 setexponent16(st0_ptr, exponent + EXTENDED_Ebias); /* Result is positive. */ 221 222 } 223