diff options
author | Linus Torvalds <torvalds@ppc970.osdl.org> | 2005-04-16 18:20:36 -0400 |
---|---|---|
committer | Linus Torvalds <torvalds@ppc970.osdl.org> | 2005-04-16 18:20:36 -0400 |
commit | 1da177e4c3f41524e886b7f1b8a0c1fc7321cac2 (patch) | |
tree | 0bba044c4ce775e45a88a51686b5d9f90697ea9d /arch/i386/math-emu/poly_tan.c |
Linux-2.6.12-rc2v2.6.12-rc2
Initial git repository build. I'm not bothering with the full history,
even though we have it. We can create a separate "historical" git
archive of that later if we want to, and in the meantime it's about
3.2GB when imported into git - space that would just make the early
git days unnecessarily complicated, when we don't have a lot of good
infrastructure for it.
Let it rip!
Diffstat (limited to 'arch/i386/math-emu/poly_tan.c')
-rw-r--r-- | arch/i386/math-emu/poly_tan.c | 222 |
1 files changed, 222 insertions, 0 deletions
diff --git a/arch/i386/math-emu/poly_tan.c b/arch/i386/math-emu/poly_tan.c new file mode 100644 index 000000000000..8df3e03b6e6f --- /dev/null +++ b/arch/i386/math-emu/poly_tan.c | |||
@@ -0,0 +1,222 @@ | |||
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 | |||
20 | |||
21 | #define HiPOWERop 3 /* odd poly, positive terms */ | ||
22 | static const unsigned long long oddplterm[HiPOWERop] = | ||
23 | { | ||
24 | 0x0000000000000000LL, | ||
25 | 0x0051a1cf08fca228LL, | ||
26 | 0x0000000071284ff7LL | ||
27 | }; | ||
28 | |||
29 | #define HiPOWERon 2 /* odd poly, negative terms */ | ||
30 | static const unsigned long long oddnegterm[HiPOWERon] = | ||
31 | { | ||
32 | 0x1291a9a184244e80LL, | ||
33 | 0x0000583245819c21LL | ||
34 | }; | ||
35 | |||
36 | #define HiPOWERep 2 /* even poly, positive terms */ | ||
37 | static const unsigned long long evenplterm[HiPOWERep] = | ||
38 | { | ||
39 | 0x0e848884b539e888LL, | ||
40 | 0x00003c7f18b887daLL | ||
41 | }; | ||
42 | |||
43 | #define HiPOWERen 2 /* even poly, negative terms */ | ||
44 | static const unsigned long long evennegterm[HiPOWERen] = | ||
45 | { | ||
46 | 0xf1f0200fd51569ccLL, | ||
47 | 0x003afb46105c4432LL | ||
48 | }; | ||
49 | |||
50 | static const unsigned long long twothirds = 0xaaaaaaaaaaaaaaabLL; | ||
51 | |||
52 | |||
53 | /*--- poly_tan() ------------------------------------------------------------+ | ||
54 | | | | ||
55 | +---------------------------------------------------------------------------*/ | ||
56 | void poly_tan(FPU_REG *st0_ptr) | ||
57 | { | ||
58 | long int exponent; | ||
59 | int invert; | ||
60 | Xsig argSq, argSqSq, accumulatoro, accumulatore, accum, | ||
61 | argSignif, fix_up; | ||
62 | unsigned long adj; | ||
63 | |||
64 | exponent = exponent(st0_ptr); | ||
65 | |||
66 | #ifdef PARANOID | ||
67 | if ( signnegative(st0_ptr) ) /* Can't hack a number < 0.0 */ | ||
68 | { arith_invalid(0); return; } /* Need a positive number */ | ||
69 | #endif /* PARANOID */ | ||
70 | |||
71 | /* Split the problem into two domains, smaller and larger than pi/4 */ | ||
72 | if ( (exponent == 0) || ((exponent == -1) && (st0_ptr->sigh > 0xc90fdaa2)) ) | ||
73 | { | ||
74 | /* The argument is greater than (approx) pi/4 */ | ||
75 | invert = 1; | ||
76 | accum.lsw = 0; | ||
77 | XSIG_LL(accum) = significand(st0_ptr); | ||
78 | |||
79 | if ( exponent == 0 ) | ||
80 | { | ||
81 | /* The argument is >= 1.0 */ | ||
82 | /* Put the binary point at the left. */ | ||
83 | XSIG_LL(accum) <<= 1; | ||
84 | } | ||
85 | /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */ | ||
86 | XSIG_LL(accum) = 0x921fb54442d18469LL - XSIG_LL(accum); | ||
87 | /* This is a special case which arises due to rounding. */ | ||
88 | if ( XSIG_LL(accum) == 0xffffffffffffffffLL ) | ||
89 | { | ||
90 | FPU_settag0(TAG_Valid); | ||
91 | significand(st0_ptr) = 0x8a51e04daabda360LL; | ||
92 | setexponent16(st0_ptr, (0x41 + EXTENDED_Ebias) | SIGN_Negative); | ||
93 | return; | ||
94 | } | ||
95 | |||
96 | argSignif.lsw = accum.lsw; | ||
97 | XSIG_LL(argSignif) = XSIG_LL(accum); | ||
98 | exponent = -1 + norm_Xsig(&argSignif); | ||
99 | } | ||
100 | else | ||
101 | { | ||
102 | invert = 0; | ||
103 | argSignif.lsw = 0; | ||
104 | XSIG_LL(accum) = XSIG_LL(argSignif) = significand(st0_ptr); | ||
105 | |||
106 | if ( exponent < -1 ) | ||
107 | { | ||
108 | /* shift the argument right by the required places */ | ||
109 | if ( FPU_shrx(&XSIG_LL(accum), -1-exponent) >= 0x80000000U ) | ||
110 | XSIG_LL(accum) ++; /* round up */ | ||
111 | } | ||
112 | } | ||
113 | |||
114 | XSIG_LL(argSq) = XSIG_LL(accum); argSq.lsw = accum.lsw; | ||
115 | mul_Xsig_Xsig(&argSq, &argSq); | ||
116 | XSIG_LL(argSqSq) = XSIG_LL(argSq); argSqSq.lsw = argSq.lsw; | ||
117 | mul_Xsig_Xsig(&argSqSq, &argSqSq); | ||
118 | |||
119 | /* Compute the negative terms for the numerator polynomial */ | ||
120 | accumulatoro.msw = accumulatoro.midw = accumulatoro.lsw = 0; | ||
121 | polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddnegterm, HiPOWERon-1); | ||
122 | mul_Xsig_Xsig(&accumulatoro, &argSq); | ||
123 | negate_Xsig(&accumulatoro); | ||
124 | /* Add the positive terms */ | ||
125 | polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddplterm, HiPOWERop-1); | ||
126 | |||
127 | |||
128 | /* Compute the positive terms for the denominator polynomial */ | ||
129 | accumulatore.msw = accumulatore.midw = accumulatore.lsw = 0; | ||
130 | polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evenplterm, HiPOWERep-1); | ||
131 | mul_Xsig_Xsig(&accumulatore, &argSq); | ||
132 | negate_Xsig(&accumulatore); | ||
133 | /* Add the negative terms */ | ||
134 | polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evennegterm, HiPOWERen-1); | ||
135 | /* Multiply by arg^2 */ | ||
136 | mul64_Xsig(&accumulatore, &XSIG_LL(argSignif)); | ||
137 | mul64_Xsig(&accumulatore, &XSIG_LL(argSignif)); | ||
138 | /* de-normalize and divide by 2 */ | ||
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 | } | ||
190 | else | ||
191 | adj = 0; | ||
192 | adj = mul_32_32(0x898cc517, adj); /* delta * tan^2 */ | ||
193 | |||
194 | fix_up.msw += adj; | ||
195 | if ( !(fix_up.msw & 0x80000000) ) /* did fix_up overflow ? */ | ||
196 | { | ||
197 | /* Yes, we need to add an msb */ | ||
198 | shr_Xsig(&fix_up, 1); | ||
199 | fix_up.msw |= 0x80000000; | ||
200 | shr_Xsig(&fix_up, 64 + exponent); | ||
201 | } | ||
202 | else | ||
203 | shr_Xsig(&fix_up, 65 + exponent); | ||
204 | |||
205 | add_two_Xsig(&accum, &fix_up, &exponent); | ||
206 | |||
207 | /* accum now contains tan(pi/2 - arg). | ||
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 | |||
222 | } | ||