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