aboutsummaryrefslogtreecommitdiffstats
path: root/arch/x86/math-emu/poly_tan.c
diff options
context:
space:
mode:
Diffstat (limited to 'arch/x86/math-emu/poly_tan.c')
-rw-r--r--arch/x86/math-emu/poly_tan.c338
1 files changed, 164 insertions, 174 deletions
diff --git a/arch/x86/math-emu/poly_tan.c b/arch/x86/math-emu/poly_tan.c
index 8df3e03b6e6f..c0d181e39229 100644
--- a/arch/x86/math-emu/poly_tan.c
+++ b/arch/x86/math-emu/poly_tan.c
@@ -17,206 +17,196 @@
17#include "control_w.h" 17#include "control_w.h"
18#include "poly.h" 18#include "poly.h"
19 19
20
21#define HiPOWERop 3 /* odd poly, positive terms */ 20#define HiPOWERop 3 /* odd poly, positive terms */
22static const unsigned long long oddplterm[HiPOWERop] = 21static const unsigned long long oddplterm[HiPOWERop] = {
23{ 22 0x0000000000000000LL,
24 0x0000000000000000LL, 23 0x0051a1cf08fca228LL,
25 0x0051a1cf08fca228LL, 24 0x0000000071284ff7LL
26 0x0000000071284ff7LL
27}; 25};
28 26
29#define HiPOWERon 2 /* odd poly, negative terms */ 27#define HiPOWERon 2 /* odd poly, negative terms */
30static const unsigned long long oddnegterm[HiPOWERon] = 28static const unsigned long long oddnegterm[HiPOWERon] = {
31{ 29 0x1291a9a184244e80LL,
32 0x1291a9a184244e80LL, 30 0x0000583245819c21LL
33 0x0000583245819c21LL
34}; 31};
35 32
36#define HiPOWERep 2 /* even poly, positive terms */ 33#define HiPOWERep 2 /* even poly, positive terms */
37static const unsigned long long evenplterm[HiPOWERep] = 34static const unsigned long long evenplterm[HiPOWERep] = {
38{ 35 0x0e848884b539e888LL,
39 0x0e848884b539e888LL, 36 0x00003c7f18b887daLL
40 0x00003c7f18b887daLL
41}; 37};
42 38
43#define HiPOWERen 2 /* even poly, negative terms */ 39#define HiPOWERen 2 /* even poly, negative terms */
44static const unsigned long long evennegterm[HiPOWERen] = 40static const unsigned long long evennegterm[HiPOWERen] = {
45{ 41 0xf1f0200fd51569ccLL,
46 0xf1f0200fd51569ccLL, 42 0x003afb46105c4432LL
47 0x003afb46105c4432LL
48}; 43};
49 44
50static const unsigned long long twothirds = 0xaaaaaaaaaaaaaaabLL; 45static const unsigned long long twothirds = 0xaaaaaaaaaaaaaaabLL;
51 46
52
53/*--- poly_tan() ------------------------------------------------------------+ 47/*--- poly_tan() ------------------------------------------------------------+
54 | | 48 | |
55 +---------------------------------------------------------------------------*/ 49 +---------------------------------------------------------------------------*/
56void poly_tan(FPU_REG *st0_ptr) 50void poly_tan(FPU_REG * st0_ptr)
57{ 51{
58 long int exponent; 52 long int exponent;
59 int invert; 53 int invert;
60 Xsig argSq, argSqSq, accumulatoro, accumulatore, accum, 54 Xsig argSq, argSqSq, accumulatoro, accumulatore, accum,
61 argSignif, fix_up; 55 argSignif, fix_up;
62 unsigned long adj; 56 unsigned long adj;
63 57
64 exponent = exponent(st0_ptr); 58 exponent = exponent(st0_ptr);
65 59
66#ifdef PARANOID 60#ifdef PARANOID
67 if ( signnegative(st0_ptr) ) /* Can't hack a number < 0.0 */ 61 if (signnegative(st0_ptr)) { /* Can't hack a number < 0.0 */
68 { arith_invalid(0); return; } /* Need a positive number */ 62 arith_invalid(0);
63 return;
64 } /* Need a positive number */
69#endif /* PARANOID */ 65#endif /* PARANOID */
70 66
71 /* Split the problem into two domains, smaller and larger than pi/4 */ 67 /* Split the problem into two domains, smaller and larger than pi/4 */
72 if ( (exponent == 0) || ((exponent == -1) && (st0_ptr->sigh > 0xc90fdaa2)) ) 68 if ((exponent == 0)
73 { 69 || ((exponent == -1) && (st0_ptr->sigh > 0xc90fdaa2))) {
74 /* The argument is greater than (approx) pi/4 */ 70 /* The argument is greater than (approx) pi/4 */
75 invert = 1; 71 invert = 1;
76 accum.lsw = 0; 72 accum.lsw = 0;
77 XSIG_LL(accum) = significand(st0_ptr); 73 XSIG_LL(accum) = significand(st0_ptr);
78 74
79 if ( exponent == 0 ) 75 if (exponent == 0) {
80 { 76 /* The argument is >= 1.0 */
81 /* The argument is >= 1.0 */ 77 /* Put the binary point at the left. */
82 /* Put the binary point at the left. */ 78 XSIG_LL(accum) <<= 1;
83 XSIG_LL(accum) <<= 1; 79 }
84 } 80 /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
85 /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */ 81 XSIG_LL(accum) = 0x921fb54442d18469LL - XSIG_LL(accum);
86 XSIG_LL(accum) = 0x921fb54442d18469LL - XSIG_LL(accum); 82 /* This is a special case which arises due to rounding. */
87 /* This is a special case which arises due to rounding. */ 83 if (XSIG_LL(accum) == 0xffffffffffffffffLL) {
88 if ( XSIG_LL(accum) == 0xffffffffffffffffLL ) 84 FPU_settag0(TAG_Valid);
89 { 85 significand(st0_ptr) = 0x8a51e04daabda360LL;
90 FPU_settag0(TAG_Valid); 86 setexponent16(st0_ptr,
91 significand(st0_ptr) = 0x8a51e04daabda360LL; 87 (0x41 + EXTENDED_Ebias) | SIGN_Negative);
92 setexponent16(st0_ptr, (0x41 + EXTENDED_Ebias) | SIGN_Negative); 88 return;
93 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 }
94 } 105 }
95 106
96 argSignif.lsw = accum.lsw; 107 XSIG_LL(argSq) = XSIG_LL(accum);
97 XSIG_LL(argSignif) = XSIG_LL(accum); 108 argSq.lsw = accum.lsw;
98 exponent = -1 + norm_Xsig(&argSignif); 109 mul_Xsig_Xsig(&argSq, &argSq);
99 } 110 XSIG_LL(argSqSq) = XSIG_LL(argSq);
100 else 111 argSqSq.lsw = argSq.lsw;
101 { 112 mul_Xsig_Xsig(&argSqSq, &argSqSq);
102 invert = 0; 113
103 argSignif.lsw = 0; 114 /* Compute the negative terms for the numerator polynomial */
104 XSIG_LL(accum) = XSIG_LL(argSignif) = significand(st0_ptr); 115 accumulatoro.msw = accumulatoro.midw = accumulatoro.lsw = 0;
105 116 polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddnegterm,
106 if ( exponent < -1 ) 117 HiPOWERon - 1);
107 { 118 mul_Xsig_Xsig(&accumulatoro, &argSq);
108 /* shift the argument right by the required places */ 119 negate_Xsig(&accumulatoro);
109 if ( FPU_shrx(&XSIG_LL(accum), -1-exponent) >= 0x80000000U ) 120 /* Add the positive terms */
110 XSIG_LL(accum) ++; /* round up */ 121 polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddplterm,
111 } 122 HiPOWERop - 1);
112 } 123
113 124 /* Compute the positive terms for the denominator polynomial */
114 XSIG_LL(argSq) = XSIG_LL(accum); argSq.lsw = accum.lsw; 125 accumulatore.msw = accumulatore.midw = accumulatore.lsw = 0;
115 mul_Xsig_Xsig(&argSq, &argSq); 126 polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evenplterm,
116 XSIG_LL(argSqSq) = XSIG_LL(argSq); argSqSq.lsw = argSq.lsw; 127 HiPOWERep - 1);
117 mul_Xsig_Xsig(&argSqSq, &argSqSq); 128 mul_Xsig_Xsig(&accumulatore, &argSq);
118 129 negate_Xsig(&accumulatore);
119 /* Compute the negative terms for the numerator polynomial */ 130 /* Add the negative terms */
120 accumulatoro.msw = accumulatoro.midw = accumulatoro.lsw = 0; 131 polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evennegterm,
121 polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddnegterm, HiPOWERon-1); 132 HiPOWERen - 1);
122 mul_Xsig_Xsig(&accumulatoro, &argSq); 133 /* Multiply by arg^2 */
123 negate_Xsig(&accumulatoro); 134 mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
124 /* Add the positive terms */ 135 mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
125 polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddplterm, HiPOWERop-1); 136 /* de-normalize and divide by 2 */
126 137 shr_Xsig(&accumulatore, -2 * (1 + exponent) + 1);
127 138 negate_Xsig(&accumulatore); /* This does 1 - accumulator */
128 /* Compute the positive terms for the denominator polynomial */ 139
129 accumulatore.msw = accumulatore.midw = accumulatore.lsw = 0; 140 /* Now find the ratio. */
130 polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evenplterm, HiPOWERep-1); 141 if (accumulatore.msw == 0) {
131 mul_Xsig_Xsig(&accumulatore, &argSq); 142 /* accumulatoro must contain 1.0 here, (actually, 0) but it
132 negate_Xsig(&accumulatore); 143 really doesn't matter what value we use because it will
133 /* Add the negative terms */ 144 have negligible effect in later calculations
134 polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evennegterm, HiPOWERen-1); 145 */
135 /* Multiply by arg^2 */ 146 XSIG_LL(accum) = 0x8000000000000000LL;
136 mul64_Xsig(&accumulatore, &XSIG_LL(argSignif)); 147 accum.lsw = 0;
137 mul64_Xsig(&accumulatore, &XSIG_LL(argSignif)); 148 } else {
138 /* de-normalize and divide by 2 */ 149 div_Xsig(&accumulatoro, &accumulatore, &accum);
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 } 150 }
190 else 151
191 adj = 0; 152 /* Multiply by 1/3 * arg^3 */
192 adj = mul_32_32(0x898cc517, adj); /* delta * tan^2 */ 153 mul64_Xsig(&accum, &XSIG_LL(argSignif));
193 154 mul64_Xsig(&accum, &XSIG_LL(argSignif));
194 fix_up.msw += adj; 155 mul64_Xsig(&accum, &XSIG_LL(argSignif));
195 if ( !(fix_up.msw & 0x80000000) ) /* did fix_up overflow ? */ 156 mul64_Xsig(&accum, &twothirds);
196 { 157 shr_Xsig(&accum, -2 * (exponent + 1));
197 /* Yes, we need to add an msb */ 158
198 shr_Xsig(&fix_up, 1); 159 /* tan(arg) = arg + accum */
199 fix_up.msw |= 0x80000000; 160 add_two_Xsig(&accum, &argSignif, &exponent);
200 shr_Xsig(&fix_up, 64 + 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;
201 } 204 }
202 else 205
203 shr_Xsig(&fix_up, 65 + exponent); 206 /* Transfer the result */
204 207 round_Xsig(&accum);
205 add_two_Xsig(&accum, &fix_up, &exponent); 208 FPU_settag0(TAG_Valid);
206 209 significand(st0_ptr) = XSIG_LL(accum);
207 /* accum now contains tan(pi/2 - arg). 210 setexponent16(st0_ptr, exponent + EXTENDED_Ebias); /* Result is positive. */
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 211
222} 212}
generated by cgit v1.2.2 (git 2.25.0) at 2025-01-08 16:41:48 -0500