File: /usr/src/linux/arch/i386/math-emu/poly_atan.c
1 /*---------------------------------------------------------------------------+
2 | poly_atan.c |
3 | |
4 | Compute the arctan of a FPU_REG, using a polynomial approximation. |
5 | |
6 | Copyright (C) 1992,1993,1994,1997 |
7 | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
8 | E-mail billm@suburbia.net |
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 "status_w.h"
18 #include "control_w.h"
19 #include "poly.h"
20
21
22 #define HIPOWERon 6 /* odd poly, negative terms */
23 static const unsigned long long oddnegterms[HIPOWERon] =
24 {
25 0x0000000000000000LL, /* Dummy (not for - 1.0) */
26 0x015328437f756467LL,
27 0x0005dda27b73dec6LL,
28 0x0000226bf2bfb91aLL,
29 0x000000ccc439c5f7LL,
30 0x0000000355438407LL
31 } ;
32
33 #define HIPOWERop 6 /* odd poly, positive terms */
34 static const unsigned long long oddplterms[HIPOWERop] =
35 {
36 /* 0xaaaaaaaaaaaaaaabLL, transferred to fixedpterm[] */
37 0x0db55a71875c9ac2LL,
38 0x0029fce2d67880b0LL,
39 0x0000dfd3908b4596LL,
40 0x00000550fd61dab4LL,
41 0x0000001c9422b3f9LL,
42 0x000000003e3301e1LL
43 };
44
45 static const unsigned long long denomterm = 0xebd9b842c5c53a0eLL;
46
47 static const Xsig fixedpterm = MK_XSIG(0xaaaaaaaa, 0xaaaaaaaa, 0xaaaaaaaa);
48
49 static const Xsig pi_signif = MK_XSIG(0xc90fdaa2, 0x2168c234, 0xc4c6628b);
50
51
52 /*--- poly_atan() -----------------------------------------------------------+
53 | |
54 +---------------------------------------------------------------------------*/
55 void poly_atan(FPU_REG *st0_ptr, u_char st0_tag,
56 FPU_REG *st1_ptr, u_char st1_tag)
57 {
58 u_char transformed, inverted,
59 sign1, sign2;
60 int exponent;
61 long int dummy_exp;
62 Xsig accumulator, Numer, Denom, accumulatore, argSignif,
63 argSq, argSqSq;
64 u_char tag;
65
66 sign1 = getsign(st0_ptr);
67 sign2 = getsign(st1_ptr);
68 if ( st0_tag == TAG_Valid )
69 {
70 exponent = exponent(st0_ptr);
71 }
72 else
73 {
74 /* This gives non-compatible stack contents... */
75 FPU_to_exp16(st0_ptr, st0_ptr);
76 exponent = exponent16(st0_ptr);
77 }
78 if ( st1_tag == TAG_Valid )
79 {
80 exponent -= exponent(st1_ptr);
81 }
82 else
83 {
84 /* This gives non-compatible stack contents... */
85 FPU_to_exp16(st1_ptr, st1_ptr);
86 exponent -= exponent16(st1_ptr);
87 }
88
89 if ( (exponent < 0) || ((exponent == 0) &&
90 ((st0_ptr->sigh < st1_ptr->sigh) ||
91 ((st0_ptr->sigh == st1_ptr->sigh) &&
92 (st0_ptr->sigl < st1_ptr->sigl))) ) )
93 {
94 inverted = 1;
95 Numer.lsw = Denom.lsw = 0;
96 XSIG_LL(Numer) = significand(st0_ptr);
97 XSIG_LL(Denom) = significand(st1_ptr);
98 }
99 else
100 {
101 inverted = 0;
102 exponent = -exponent;
103 Numer.lsw = Denom.lsw = 0;
104 XSIG_LL(Numer) = significand(st1_ptr);
105 XSIG_LL(Denom) = significand(st0_ptr);
106 }
107 div_Xsig(&Numer, &Denom, &argSignif);
108 exponent += norm_Xsig(&argSignif);
109
110 if ( (exponent >= -1)
111 || ((exponent == -2) && (argSignif.msw > 0xd413ccd0)) )
112 {
113 /* The argument is greater than sqrt(2)-1 (=0.414213562...) */
114 /* Convert the argument by an identity for atan */
115 transformed = 1;
116
117 if ( exponent >= 0 )
118 {
119 #ifdef PARANOID
120 if ( !( (exponent == 0) &&
121 (argSignif.lsw == 0) && (argSignif.midw == 0) &&
122 (argSignif.msw == 0x80000000) ) )
123 {
124 EXCEPTION(EX_INTERNAL|0x104); /* There must be a logic error */
125 return;
126 }
127 #endif /* PARANOID */
128 argSignif.msw = 0; /* Make the transformed arg -> 0.0 */
129 }
130 else
131 {
132 Numer.lsw = Denom.lsw = argSignif.lsw;
133 XSIG_LL(Numer) = XSIG_LL(Denom) = XSIG_LL(argSignif);
134
135 if ( exponent < -1 )
136 shr_Xsig(&Numer, -1-exponent);
137 negate_Xsig(&Numer);
138
139 shr_Xsig(&Denom, -exponent);
140 Denom.msw |= 0x80000000;
141
142 div_Xsig(&Numer, &Denom, &argSignif);
143
144 exponent = -1 + norm_Xsig(&argSignif);
145 }
146 }
147 else
148 {
149 transformed = 0;
150 }
151
152 argSq.lsw = argSignif.lsw; argSq.midw = argSignif.midw;
153 argSq.msw = argSignif.msw;
154 mul_Xsig_Xsig(&argSq, &argSq);
155
156 argSqSq.lsw = argSq.lsw; argSqSq.midw = argSq.midw; argSqSq.msw = argSq.msw;
157 mul_Xsig_Xsig(&argSqSq, &argSqSq);
158
159 accumulatore.lsw = argSq.lsw;
160 XSIG_LL(accumulatore) = XSIG_LL(argSq);
161
162 shr_Xsig(&argSq, 2*(-1-exponent-1));
163 shr_Xsig(&argSqSq, 4*(-1-exponent-1));
164
165 /* Now have argSq etc with binary point at the left
166 .1xxxxxxxx */
167
168 /* Do the basic fixed point polynomial evaluation */
169 accumulator.msw = accumulator.midw = accumulator.lsw = 0;
170 polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq),
171 oddplterms, HIPOWERop-1);
172 mul64_Xsig(&accumulator, &XSIG_LL(argSq));
173 negate_Xsig(&accumulator);
174 polynomial_Xsig(&accumulator, &XSIG_LL(argSqSq), oddnegterms, HIPOWERon-1);
175 negate_Xsig(&accumulator);
176 add_two_Xsig(&accumulator, &fixedpterm, &dummy_exp);
177
178 mul64_Xsig(&accumulatore, &denomterm);
179 shr_Xsig(&accumulatore, 1 + 2*(-1-exponent));
180 accumulatore.msw |= 0x80000000;
181
182 div_Xsig(&accumulator, &accumulatore, &accumulator);
183
184 mul_Xsig_Xsig(&accumulator, &argSignif);
185 mul_Xsig_Xsig(&accumulator, &argSq);
186
187 shr_Xsig(&accumulator, 3);
188 negate_Xsig(&accumulator);
189 add_Xsig_Xsig(&accumulator, &argSignif);
190
191 if ( transformed )
192 {
193 /* compute pi/4 - accumulator */
194 shr_Xsig(&accumulator, -1-exponent);
195 negate_Xsig(&accumulator);
196 add_Xsig_Xsig(&accumulator, &pi_signif);
197 exponent = -1;
198 }
199
200 if ( inverted )
201 {
202 /* compute pi/2 - accumulator */
203 shr_Xsig(&accumulator, -exponent);
204 negate_Xsig(&accumulator);
205 add_Xsig_Xsig(&accumulator, &pi_signif);
206 exponent = 0;
207 }
208
209 if ( sign1 )
210 {
211 /* compute pi - accumulator */
212 shr_Xsig(&accumulator, 1 - exponent);
213 negate_Xsig(&accumulator);
214 add_Xsig_Xsig(&accumulator, &pi_signif);
215 exponent = 1;
216 }
217
218 exponent += round_Xsig(&accumulator);
219
220 significand(st1_ptr) = XSIG_LL(accumulator);
221 setexponent16(st1_ptr, exponent);
222
223 tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign2);
224 FPU_settagi(1, tag);
225
226 set_precision_flag_up(); /* We do not really know if up or down,
227 use this as the default. */
228
229 }
230