← Back to branch summary

~ubuntu-branches/ubuntu/precise/linux-lowlatency/precise

« back to all changes in this revision

Viewing changes to arch/x86/math-emu/poly_sin.c

  • Committer: Package Import Robot
  • Author(s): Alessio Igor Bogani
  • Date: 2011-10-26 11:13:05 UTC
  • Revision ID: package-import@ubuntu.com-20111026111305-tz023xykf0i6eosh
Tags: upstream-3.2.0
Import upstream version 3.2.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
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
}