File: /usr/src/linux/arch/i386/math-emu/poly_tan.c
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 }
223