aboutsummaryrefslogtreecommitdiffstats
path: root/arch/x86/math-emu/poly_tan.c
diff options
context:
space:
mode:
authorIngo Molnar <mingo@elte.hu>2008-01-30 07:30:11 -0500
committerIngo Molnar <mingo@elte.hu>2008-01-30 07:30:11 -0500
commit3d0d14f983b55a570b976976284df4c434af3223 (patch)
tree864f11c0ce5ee1e15acdd196018b79d0d0e2685d /arch/x86/math-emu/poly_tan.c
parenta4ec1effce83796209a0258602b0cf50026d86f2 (diff)
x86: lindent arch/i386/math-emu
lindent these files: errors lines of code errors/KLOC arch/x86/math-emu/ 2236 9424 237.2 arch/x86/math-emu/ 128 8706 14.7 no other changes. No code changed: text data bss dec hex filename 5589802 612739 3833856 10036397 9924ad vmlinux.before 5589802 612739 3833856 10036397 9924ad vmlinux.after the intent of this patch is to ease the automated tracking of kernel code quality - it's just much easier for us to maintain it if every file in arch/x86 is supposed to be clean. NOTE: it is a known problem of lindent that it causes some style damage of its own, but it's a safe tool (well, except for the gcc array range initializers extension), so we did the bulk of the changes via lindent, and did the manual fixups in a followup patch. the resulting math-emu code has been tested by Thomas Gleixner on a real 386 DX CPU as well, and it works fine. Signed-off-by: Ingo Molnar <mingo@elte.hu> Signed-off-by: Thomas Gleixner <tglx@linutronix.de>
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}