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