Merge branches 'release' and 'gpe-ack' into release
[linux-2.6] / arch / x86 / math-emu / poly_sin.c
1 /*---------------------------------------------------------------------------+
2  |  poly_sin.c                                                               |
3  |                                                                           |
4  |  Computation of an approximation of the sin function and the cosine       |
5  |  function by a polynomial.                                                |
6  |                                                                           |
7  | Copyright (C) 1992,1993,1994,1997,1999                                    |
8  |                  W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
9  |                  E-mail   billm@melbpc.org.au                             |
10  |                                                                           |
11  |                                                                           |
12  +---------------------------------------------------------------------------*/
13
14 #include "exception.h"
15 #include "reg_constant.h"
16 #include "fpu_emu.h"
17 #include "fpu_system.h"
18 #include "control_w.h"
19 #include "poly.h"
20
21 #define N_COEFF_P       4
22 #define N_COEFF_N       4
23
24 static const unsigned long long pos_terms_l[N_COEFF_P] = {
25         0xaaaaaaaaaaaaaaabLL,
26         0x00d00d00d00cf906LL,
27         0x000006b99159a8bbLL,
28         0x000000000d7392e6LL
29 };
30
31 static const unsigned long long neg_terms_l[N_COEFF_N] = {
32         0x2222222222222167LL,
33         0x0002e3bc74aab624LL,
34         0x0000000b09229062LL,
35         0x00000000000c7973LL
36 };
37
38 #define N_COEFF_PH      4
39 #define N_COEFF_NH      4
40 static const unsigned long long pos_terms_h[N_COEFF_PH] = {
41         0x0000000000000000LL,
42         0x05b05b05b05b0406LL,
43         0x000049f93edd91a9LL,
44         0x00000000c9c9ed62LL
45 };
46
47 static const unsigned long long neg_terms_h[N_COEFF_NH] = {
48         0xaaaaaaaaaaaaaa98LL,
49         0x001a01a01a019064LL,
50         0x0000008f76c68a77LL,
51         0x0000000000d58f5eLL
52 };
53
54 /*--- poly_sine() -----------------------------------------------------------+
55  |                                                                           |
56  +---------------------------------------------------------------------------*/
57 void poly_sine(FPU_REG *st0_ptr)
58 {
59         int exponent, echange;
60         Xsig accumulator, argSqrd, argTo4;
61         unsigned long fix_up, adj;
62         unsigned long long fixed_arg;
63         FPU_REG result;
64
65         exponent = exponent(st0_ptr);
66
67         accumulator.lsw = accumulator.midw = accumulator.msw = 0;
68
69         /* Split into two ranges, for arguments below and above 1.0 */
70         /* The boundary between upper and lower is approx 0.88309101259 */
71         if ((exponent < -1)
72             || ((exponent == -1) && (st0_ptr->sigh <= 0xe21240aa))) {
73                 /* The argument is <= 0.88309101259 */
74
75                 argSqrd.msw = st0_ptr->sigh;
76                 argSqrd.midw = st0_ptr->sigl;
77                 argSqrd.lsw = 0;
78                 mul64_Xsig(&argSqrd, &significand(st0_ptr));
79                 shr_Xsig(&argSqrd, 2 * (-1 - exponent));
80                 argTo4.msw = argSqrd.msw;
81                 argTo4.midw = argSqrd.midw;
82                 argTo4.lsw = argSqrd.lsw;
83                 mul_Xsig_Xsig(&argTo4, &argTo4);
84
85                 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_l,
86                                 N_COEFF_N - 1);
87                 mul_Xsig_Xsig(&accumulator, &argSqrd);
88                 negate_Xsig(&accumulator);
89
90                 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_l,
91                                 N_COEFF_P - 1);
92
93                 shr_Xsig(&accumulator, 2);      /* Divide by four */
94                 accumulator.msw |= 0x80000000;  /* Add 1.0 */
95
96                 mul64_Xsig(&accumulator, &significand(st0_ptr));
97                 mul64_Xsig(&accumulator, &significand(st0_ptr));
98                 mul64_Xsig(&accumulator, &significand(st0_ptr));
99
100                 /* Divide by four, FPU_REG compatible, etc */
101                 exponent = 3 * exponent;
102
103                 /* The minimum exponent difference is 3 */
104                 shr_Xsig(&accumulator, exponent(st0_ptr) - exponent);
105
106                 negate_Xsig(&accumulator);
107                 XSIG_LL(accumulator) += significand(st0_ptr);
108
109                 echange = round_Xsig(&accumulator);
110
111                 setexponentpos(&result, exponent(st0_ptr) + echange);
112         } else {
113                 /* The argument is > 0.88309101259 */
114                 /* We use sin(st(0)) = cos(pi/2-st(0)) */
115
116                 fixed_arg = significand(st0_ptr);
117
118                 if (exponent == 0) {
119                         /* The argument is >= 1.0 */
120
121                         /* Put the binary point at the left. */
122                         fixed_arg <<= 1;
123                 }
124                 /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
125                 fixed_arg = 0x921fb54442d18469LL - fixed_arg;
126                 /* There is a special case which arises due to rounding, to fix here. */
127                 if (fixed_arg == 0xffffffffffffffffLL)
128                         fixed_arg = 0;
129
130                 XSIG_LL(argSqrd) = fixed_arg;
131                 argSqrd.lsw = 0;
132                 mul64_Xsig(&argSqrd, &fixed_arg);
133
134                 XSIG_LL(argTo4) = XSIG_LL(argSqrd);
135                 argTo4.lsw = argSqrd.lsw;
136                 mul_Xsig_Xsig(&argTo4, &argTo4);
137
138                 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_h,
139                                 N_COEFF_NH - 1);
140                 mul_Xsig_Xsig(&accumulator, &argSqrd);
141                 negate_Xsig(&accumulator);
142
143                 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_h,
144                                 N_COEFF_PH - 1);
145                 negate_Xsig(&accumulator);
146
147                 mul64_Xsig(&accumulator, &fixed_arg);
148                 mul64_Xsig(&accumulator, &fixed_arg);
149
150                 shr_Xsig(&accumulator, 3);
151                 negate_Xsig(&accumulator);
152
153                 add_Xsig_Xsig(&accumulator, &argSqrd);
154
155                 shr_Xsig(&accumulator, 1);
156
157                 accumulator.lsw |= 1;   /* A zero accumulator here would cause problems */
158                 negate_Xsig(&accumulator);
159
160                 /* The basic computation is complete. Now fix the answer to
161                    compensate for the error due to the approximation used for
162                    pi/2
163                  */
164
165                 /* This has an exponent of -65 */
166                 fix_up = 0x898cc517;
167                 /* The fix-up needs to be improved for larger args */
168                 if (argSqrd.msw & 0xffc00000) {
169                         /* Get about 32 bit precision in these: */
170                         fix_up -= mul_32_32(0x898cc517, argSqrd.msw) / 6;
171                 }
172                 fix_up = mul_32_32(fix_up, LL_MSW(fixed_arg));
173
174                 adj = accumulator.lsw;  /* temp save */
175                 accumulator.lsw -= fix_up;
176                 if (accumulator.lsw > adj)
177                         XSIG_LL(accumulator)--;
178
179                 echange = round_Xsig(&accumulator);
180
181                 setexponentpos(&result, echange - 1);
182         }
183
184         significand(&result) = XSIG_LL(accumulator);
185         setsign(&result, getsign(st0_ptr));
186         FPU_copy_to_reg0(&result, TAG_Valid);
187
188 #ifdef PARANOID
189         if ((exponent(&result) >= 0)
190             && (significand(&result) > 0x8000000000000000LL)) {
191                 EXCEPTION(EX_INTERNAL | 0x150);
192         }
193 #endif /* PARANOID */
194
195 }
196
197 /*--- poly_cos() ------------------------------------------------------------+
198  |                                                                           |
199  +---------------------------------------------------------------------------*/
200 void poly_cos(FPU_REG *st0_ptr)
201 {
202         FPU_REG result;
203         long int exponent, exp2, echange;
204         Xsig accumulator, argSqrd, fix_up, argTo4;
205         unsigned long long fixed_arg;
206
207 #ifdef PARANOID
208         if ((exponent(st0_ptr) > 0)
209             || ((exponent(st0_ptr) == 0)
210                 && (significand(st0_ptr) > 0xc90fdaa22168c234LL))) {
211                 EXCEPTION(EX_Invalid);
212                 FPU_copy_to_reg0(&CONST_QNaN, TAG_Special);
213                 return;
214         }
215 #endif /* PARANOID */
216
217         exponent = exponent(st0_ptr);
218
219         accumulator.lsw = accumulator.midw = accumulator.msw = 0;
220
221         if ((exponent < -1)
222             || ((exponent == -1) && (st0_ptr->sigh <= 0xb00d6f54))) {
223                 /* arg is < 0.687705 */
224
225                 argSqrd.msw = st0_ptr->sigh;
226                 argSqrd.midw = st0_ptr->sigl;
227                 argSqrd.lsw = 0;
228                 mul64_Xsig(&argSqrd, &significand(st0_ptr));
229
230                 if (exponent < -1) {
231                         /* shift the argument right by the required places */
232                         shr_Xsig(&argSqrd, 2 * (-1 - exponent));
233                 }
234
235                 argTo4.msw = argSqrd.msw;
236                 argTo4.midw = argSqrd.midw;
237                 argTo4.lsw = argSqrd.lsw;
238                 mul_Xsig_Xsig(&argTo4, &argTo4);
239
240                 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_h,
241                                 N_COEFF_NH - 1);
242                 mul_Xsig_Xsig(&accumulator, &argSqrd);
243                 negate_Xsig(&accumulator);
244
245                 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_h,
246                                 N_COEFF_PH - 1);
247                 negate_Xsig(&accumulator);
248
249                 mul64_Xsig(&accumulator, &significand(st0_ptr));
250                 mul64_Xsig(&accumulator, &significand(st0_ptr));
251                 shr_Xsig(&accumulator, -2 * (1 + exponent));
252
253                 shr_Xsig(&accumulator, 3);
254                 negate_Xsig(&accumulator);
255
256                 add_Xsig_Xsig(&accumulator, &argSqrd);
257
258                 shr_Xsig(&accumulator, 1);
259
260                 /* It doesn't matter if accumulator is all zero here, the
261                    following code will work ok */
262                 negate_Xsig(&accumulator);
263
264                 if (accumulator.lsw & 0x80000000)
265                         XSIG_LL(accumulator)++;
266                 if (accumulator.msw == 0) {
267                         /* The result is 1.0 */
268                         FPU_copy_to_reg0(&CONST_1, TAG_Valid);
269                         return;
270                 } else {
271                         significand(&result) = XSIG_LL(accumulator);
272
273                         /* will be a valid positive nr with expon = -1 */
274                         setexponentpos(&result, -1);
275                 }
276         } else {
277                 fixed_arg = significand(st0_ptr);
278
279                 if (exponent == 0) {
280                         /* The argument is >= 1.0 */
281
282                         /* Put the binary point at the left. */
283                         fixed_arg <<= 1;
284                 }
285                 /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
286                 fixed_arg = 0x921fb54442d18469LL - fixed_arg;
287                 /* There is a special case which arises due to rounding, to fix here. */
288                 if (fixed_arg == 0xffffffffffffffffLL)
289                         fixed_arg = 0;
290
291                 exponent = -1;
292                 exp2 = -1;
293
294                 /* A shift is needed here only for a narrow range of arguments,
295                    i.e. for fixed_arg approx 2^-32, but we pick up more... */
296                 if (!(LL_MSW(fixed_arg) & 0xffff0000)) {
297                         fixed_arg <<= 16;
298                         exponent -= 16;
299                         exp2 -= 16;
300                 }
301
302                 XSIG_LL(argSqrd) = fixed_arg;
303                 argSqrd.lsw = 0;
304                 mul64_Xsig(&argSqrd, &fixed_arg);
305
306                 if (exponent < -1) {
307                         /* shift the argument right by the required places */
308                         shr_Xsig(&argSqrd, 2 * (-1 - exponent));
309                 }
310
311                 argTo4.msw = argSqrd.msw;
312                 argTo4.midw = argSqrd.midw;
313                 argTo4.lsw = argSqrd.lsw;
314                 mul_Xsig_Xsig(&argTo4, &argTo4);
315
316                 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_l,
317                                 N_COEFF_N - 1);
318                 mul_Xsig_Xsig(&accumulator, &argSqrd);
319                 negate_Xsig(&accumulator);
320
321                 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_l,
322                                 N_COEFF_P - 1);
323
324                 shr_Xsig(&accumulator, 2);      /* Divide by four */
325                 accumulator.msw |= 0x80000000;  /* Add 1.0 */
326
327                 mul64_Xsig(&accumulator, &fixed_arg);
328                 mul64_Xsig(&accumulator, &fixed_arg);
329                 mul64_Xsig(&accumulator, &fixed_arg);
330
331                 /* Divide by four, FPU_REG compatible, etc */
332                 exponent = 3 * exponent;
333
334                 /* The minimum exponent difference is 3 */
335                 shr_Xsig(&accumulator, exp2 - exponent);
336
337                 negate_Xsig(&accumulator);
338                 XSIG_LL(accumulator) += fixed_arg;
339
340                 /* The basic computation is complete. Now fix the answer to
341                    compensate for the error due to the approximation used for
342                    pi/2
343                  */
344
345                 /* This has an exponent of -65 */
346                 XSIG_LL(fix_up) = 0x898cc51701b839a2ll;
347                 fix_up.lsw = 0;
348
349                 /* The fix-up needs to be improved for larger args */
350                 if (argSqrd.msw & 0xffc00000) {
351                         /* Get about 32 bit precision in these: */
352                         fix_up.msw -= mul_32_32(0x898cc517, argSqrd.msw) / 2;
353                         fix_up.msw += mul_32_32(0x898cc517, argTo4.msw) / 24;
354                 }
355
356                 exp2 += norm_Xsig(&accumulator);
357                 shr_Xsig(&accumulator, 1);      /* Prevent overflow */
358                 exp2++;
359                 shr_Xsig(&fix_up, 65 + exp2);
360
361                 add_Xsig_Xsig(&accumulator, &fix_up);
362
363                 echange = round_Xsig(&accumulator);
364
365                 setexponentpos(&result, exp2 + echange);
366                 significand(&result) = XSIG_LL(accumulator);
367         }
368
369         FPU_copy_to_reg0(&result, TAG_Valid);
370
371 #ifdef PARANOID
372         if ((exponent(&result) >= 0)
373             && (significand(&result) > 0x8000000000000000LL)) {
374                 EXCEPTION(EX_INTERNAL | 0x151);
375         }
376 #endif /* PARANOID */
377
378 }