1 | // SPDX-License-Identifier: GPL-2.0 |
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 | |
21 | #define HiPOWERop 3 /* odd poly, positive terms */ |
22 | static const unsigned long long oddplterm[HiPOWERop] = { |
23 | 0x0000000000000000LL, |
24 | 0x0051a1cf08fca228LL, |
25 | 0x0000000071284ff7LL |
26 | }; |
27 | |
28 | #define HiPOWERon 2 /* odd poly, negative terms */ |
29 | static const unsigned long long oddnegterm[HiPOWERon] = { |
30 | 0x1291a9a184244e80LL, |
31 | 0x0000583245819c21LL |
32 | }; |
33 | |
34 | #define HiPOWERep 2 /* even poly, positive terms */ |
35 | static const unsigned long long evenplterm[HiPOWERep] = { |
36 | 0x0e848884b539e888LL, |
37 | 0x00003c7f18b887daLL |
38 | }; |
39 | |
40 | #define HiPOWERen 2 /* even poly, negative terms */ |
41 | static const unsigned long long evennegterm[HiPOWERen] = { |
42 | 0xf1f0200fd51569ccLL, |
43 | 0x003afb46105c4432LL |
44 | }; |
45 | |
46 | static const unsigned long long twothirds = 0xaaaaaaaaaaaaaaabLL; |
47 | |
48 | /*--- poly_tan() ------------------------------------------------------------+ |
49 | | | |
50 | +---------------------------------------------------------------------------*/ |
51 | void poly_tan(FPU_REG *st0_ptr) |
52 | { |
53 | long int exponent; |
54 | int invert; |
55 | Xsig argSq, argSqSq, accumulatoro, accumulatore, accum, |
56 | argSignif, fix_up; |
57 | unsigned long adj; |
58 | |
59 | exponent = exponent(st0_ptr); |
60 | |
61 | #ifdef PARANOID |
62 | if (signnegative(st0_ptr)) { /* Can't hack a number < 0.0 */ |
63 | arith_invalid(0); |
64 | return; |
65 | } /* Need a positive number */ |
66 | #endif /* PARANOID */ |
67 | |
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(l: &XSIG_LL(accum), x: -1 - exponent) >= |
103 | 0x80000000U) |
104 | XSIG_LL(accum)++; /* round up */ |
105 | } |
106 | } |
107 | |
108 | XSIG_LL(argSq) = XSIG_LL(accum); |
109 | argSq.lsw = accum.lsw; |
110 | mul_Xsig_Xsig(dest: &argSq, mult: &argSq); |
111 | XSIG_LL(argSqSq) = XSIG_LL(argSq); |
112 | argSqSq.lsw = argSq.lsw; |
113 | mul_Xsig_Xsig(dest: &argSqSq, mult: &argSqSq); |
114 | |
115 | /* Compute the negative terms for the numerator polynomial */ |
116 | accumulatoro.msw = accumulatoro.midw = accumulatoro.lsw = 0; |
117 | polynomial_Xsig(&accumulatoro, x: &XSIG_LL(argSqSq), terms: oddnegterm, |
118 | HiPOWERon - 1); |
119 | mul_Xsig_Xsig(dest: &accumulatoro, mult: &argSq); |
120 | negate_Xsig(x: &accumulatoro); |
121 | /* Add the positive terms */ |
122 | polynomial_Xsig(&accumulatoro, x: &XSIG_LL(argSqSq), terms: 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, x: &XSIG_LL(argSqSq), terms: evenplterm, |
128 | HiPOWERep - 1); |
129 | mul_Xsig_Xsig(dest: &accumulatore, mult: &argSq); |
130 | negate_Xsig(x: &accumulatore); |
131 | /* Add the negative terms */ |
132 | polynomial_Xsig(&accumulatore, x: &XSIG_LL(argSqSq), terms: evennegterm, |
133 | HiPOWERen - 1); |
134 | /* Multiply by arg^2 */ |
135 | mul64_Xsig(&accumulatore, mult: &XSIG_LL(argSignif)); |
136 | mul64_Xsig(&accumulatore, mult: &XSIG_LL(argSignif)); |
137 | /* de-normalize and divide by 2 */ |
138 | shr_Xsig(&accumulatore, n: -2 * (1 + exponent) + 1); |
139 | negate_Xsig(x: &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(x1: &accumulatoro, x2: &accumulatore, dest: &accum); |
151 | } |
152 | |
153 | /* Multiply by 1/3 * arg^3 */ |
154 | mul64_Xsig(&accum, mult: &XSIG_LL(argSignif)); |
155 | mul64_Xsig(&accum, mult: &XSIG_LL(argSignif)); |
156 | mul64_Xsig(&accum, mult: &XSIG_LL(argSignif)); |
157 | mul64_Xsig(&accum, mult: &twothirds); |
158 | shr_Xsig(&accum, n: -2 * (exponent + 1)); |
159 | |
160 | /* tan(arg) = arg + accum */ |
161 | add_two_Xsig(dest: &accum, x2: &argSignif, exp: &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(arg1: adj, arg2: adj); /* tan^2 */ |
183 | } else |
184 | adj = 0; |
185 | adj = mul_32_32(arg1: 0x898cc517, arg2: 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, n: 1); |
191 | fix_up.msw |= 0x80000000; |
192 | shr_Xsig(&fix_up, n: 64 + exponent); |
193 | } else |
194 | shr_Xsig(&fix_up, n: 65 + exponent); |
195 | |
196 | add_two_Xsig(dest: &accum, x2: &fix_up, exp: &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(x1: &accumulatoro, x2: &accum, dest: &accum); |
204 | exponent = -exponent - 1; |
205 | } |
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. */ |
212 | |
213 | } |
214 | |