![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to fsincos.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Previous NeXT emulator] / previous / src / softfloat|
Return to fsincos.c CVS log | Up to [Previous NeXT emulator] / previous / src / softfloat |
1.1 root 1: /*============================================================================
2: This source file is an extension to the SoftFloat IEC/IEEE Floating-point
3: Arithmetic Package, Release 2b, written for Bochs (x86 achitecture simulator)
4: floating point emulation.
5:
6: THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE. Although reasonable effort has
7: been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT TIMES
8: RESULT IN INCORRECT BEHAVIOR. USE OF THIS SOFTWARE IS RESTRICTED TO PERSONS
9: AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ALL LOSSES,
10: COSTS, OR OTHER PROBLEMS THEY INCUR DUE TO THE SOFTWARE, AND WHO FURTHERMORE
11: EFFECTIVELY INDEMNIFY JOHN HAUSER AND THE INTERNATIONAL COMPUTER SCIENCE
12: INSTITUTE (possibly via similar legal warning) AGAINST ALL LOSSES, COSTS, OR
13: OTHER PROBLEMS INCURRED BY THEIR CUSTOMERS AND CLIENTS DUE TO THE SOFTWARE.
14:
15: Derivative works are acceptable, even for commercial purposes, so long as
16: (1) the source code for the derivative work includes prominent notice that
17: the work is derivative, and (2) the source code includes prominent notice with
18: these four paragraphs for those parts of this code that are retained.
19: =============================================================================*/
20:
21: /*============================================================================
22: * Written for Bochs (x86 achitecture simulator) by
23: * Stanislav Shwartsman [sshwarts at sourceforge net]
24: * ==========================================================================*/
25:
26: #define FLOAT128
27:
28: #define USE_estimateDiv128To64
29: #include "mamesf.h"
30: #include "softfloat.h"
31: #include "fpu_constant.h"
32:
33: static const floatx80 floatx80_one = { 0x3fff, 0x8000000000000000U };
34: static const floatx80 floatx80_default_nan = { floatx80_default_nan_high, floatx80_default_nan_low };
35:
36: #define packFloat2x128m(zHi, zLo) {(zHi), (zLo)}
37: #define PACK_FLOAT_128(hi,lo) packFloat2x128m(LIT64(hi),LIT64(lo))
38:
39: #define EXP_BIAS 0x3FFF
40:
41:
42: /* reduce trigonometric function argument using 128-bit precision
43: M_PI approximation */
44: static uint64_t argument_reduction_kernel(uint64_t aSig0, int Exp, uint64_t *zSig0, uint64_t *zSig1)
45: {
46: uint64_t term0, term1, term2;
47: uint64_t aSig1 = 0;
48:
49: shortShift128Left(aSig1, aSig0, Exp, &aSig1, &aSig0);
50: uint64_t q = estimateDiv128To64(aSig1, aSig0, FLOAT_PI_HI);
51: mul128By64To192(FLOAT_PI_HI, FLOAT_PI_LO, q, &term0, &term1, &term2);
52: sub128(aSig1, aSig0, term0, term1, zSig1, zSig0);
53: while ((int64_t)(*zSig1) < 0) {
54: --q;
55: add192(*zSig1, *zSig0, term2, 0, FLOAT_PI_HI, FLOAT_PI_LO, zSig1, zSig0, &term2);
56: }
57: *zSig1 = term2;
58: return q;
59: }
60:
61: static int reduce_trig_arg(int expDiff, int *zSign, uint64_t *aSig0, uint64_t *aSig1)
62: {
63: uint64_t term0, term1, q = 0;
64:
65: if (expDiff < 0) {
66: shift128Right(*aSig0, 0, 1, aSig0, aSig1);
67: expDiff = 0;
68: }
69: if (expDiff > 0) {
70: q = argument_reduction_kernel(*aSig0, expDiff, aSig0, aSig1);
71: }
72: else {
73: if (FLOAT_PI_HI <= *aSig0) {
74: *aSig0 -= FLOAT_PI_HI;
75: q = 1;
76: }
77: }
78:
79: shift128Right(FLOAT_PI_HI, FLOAT_PI_LO, 1, &term0, &term1);
80: if (! lt128(*aSig0, *aSig1, term0, term1))
81: {
82: int lt = lt128(term0, term1, *aSig0, *aSig1);
83: int eq = eq128(*aSig0, *aSig1, term0, term1);
84:
85: if ((eq && (q & 1)) || lt) {
86: *zSign = !(*zSign);
87: ++q;
88: }
89: if (lt) sub128(FLOAT_PI_HI, FLOAT_PI_LO, *aSig0, *aSig1, aSig0, aSig1);
90: }
91:
92: return (int)(q & 3);
93: }
94:
95: #define SIN_ARR_SIZE 11
96: #define COS_ARR_SIZE 11
97:
98: static float128 sin_arr[SIN_ARR_SIZE] =
99: {
100: PACK_FLOAT_128(0x3fff000000000000, 0x0000000000000000), /* 1 */
101: PACK_FLOAT_128(0xbffc555555555555, 0x5555555555555555), /* 3 */
102: PACK_FLOAT_128(0x3ff8111111111111, 0x1111111111111111), /* 5 */
103: PACK_FLOAT_128(0xbff2a01a01a01a01, 0xa01a01a01a01a01a), /* 7 */
104: PACK_FLOAT_128(0x3fec71de3a556c73, 0x38faac1c88e50017), /* 9 */
105: PACK_FLOAT_128(0xbfe5ae64567f544e, 0x38fe747e4b837dc7), /* 11 */
106: PACK_FLOAT_128(0x3fde6124613a86d0, 0x97ca38331d23af68), /* 13 */
107: PACK_FLOAT_128(0xbfd6ae7f3e733b81, 0xf11d8656b0ee8cb0), /* 15 */
108: PACK_FLOAT_128(0x3fce952c77030ad4, 0xa6b2605197771b00), /* 17 */
109: PACK_FLOAT_128(0xbfc62f49b4681415, 0x724ca1ec3b7b9675), /* 19 */
110: PACK_FLOAT_128(0x3fbd71b8ef6dcf57, 0x18bef146fcee6e45) /* 21 */
111: };
112:
113: static float128 cos_arr[COS_ARR_SIZE] =
114: {
115: PACK_FLOAT_128(0x3fff000000000000, 0x0000000000000000), /* 0 */
116: PACK_FLOAT_128(0xbffe000000000000, 0x0000000000000000), /* 2 */
117: PACK_FLOAT_128(0x3ffa555555555555, 0x5555555555555555), /* 4 */
118: PACK_FLOAT_128(0xbff56c16c16c16c1, 0x6c16c16c16c16c17), /* 6 */
119: PACK_FLOAT_128(0x3fefa01a01a01a01, 0xa01a01a01a01a01a), /* 8 */
120: PACK_FLOAT_128(0xbfe927e4fb7789f5, 0xc72ef016d3ea6679), /* 10 */
121: PACK_FLOAT_128(0x3fe21eed8eff8d89, 0x7b544da987acfe85), /* 12 */
122: PACK_FLOAT_128(0xbfda93974a8c07c9, 0xd20badf145dfa3e5), /* 14 */
123: PACK_FLOAT_128(0x3fd2ae7f3e733b81, 0xf11d8656b0ee8cb0), /* 16 */
124: PACK_FLOAT_128(0xbfca6827863b97d9, 0x77bb004886a2c2ab), /* 18 */
125: PACK_FLOAT_128(0x3fc1e542ba402022, 0x507a9cad2bf8f0bb) /* 20 */
126: };
127:
128: extern float128 OddPoly (float128 x, float128 *arr, unsigned n);
129:
130: /* 0 <= x <= pi/4 */
131: INLINE float128 poly_sin(float128 x)
132: {
133: // 3 5 7 9 11 13 15
134: // x x x x x x x
135: // sin (x) ~ x - --- + --- - --- + --- - ---- + ---- - ---- =
136: // 3! 5! 7! 9! 11! 13! 15!
137: //
138: // 2 4 6 8 10 12 14
139: // x x x x x x x
140: // = x * [ 1 - --- + --- - --- + --- - ---- + ---- - ---- ] =
141: // 3! 5! 7! 9! 11! 13! 15!
142: //
143: // 3 3
144: // -- 4k -- 4k+2
145: // p(x) = > C * x > 0 q(x) = > C * x < 0
146: // -- 2k -- 2k+1
147: // k=0 k=0
148: //
149: // 2
150: // sin(x) ~ x * [ p(x) + x * q(x) ]
151: //
152:
153: return OddPoly(x, sin_arr, SIN_ARR_SIZE);
154: }
155:
156: extern float128 EvenPoly(float128 x, float128 *arr, unsigned n);
157:
158: /* 0 <= x <= pi/4 */
159: INLINE float128 poly_cos(float128 x)
160: {
161: // 2 4 6 8 10 12 14
162: // x x x x x x x
163: // cos (x) ~ 1 - --- + --- - --- + --- - ---- + ---- - ----
164: // 2! 4! 6! 8! 10! 12! 14!
165: //
166: // 3 3
167: // -- 4k -- 4k+2
168: // p(x) = > C * x > 0 q(x) = > C * x < 0
169: // -- 2k -- 2k+1
170: // k=0 k=0
171: //
172: // 2
173: // cos(x) ~ [ p(x) + x * q(x) ]
174: //
175:
176: return EvenPoly(x, cos_arr, COS_ARR_SIZE);
177: }
178:
179: INLINE void sincos_invalid(floatx80 *sin_a, floatx80 *cos_a, floatx80 a)
180: {
181: if (sin_a) *sin_a = a;
182: if (cos_a) *cos_a = a;
183: }
184:
185: INLINE void sincos_tiny_argument(floatx80 *sin_a, floatx80 *cos_a, floatx80 a)
186: {
187: if (sin_a) *sin_a = a;
188: if (cos_a) *cos_a = floatx80_one;
189: }
190:
191: static floatx80 sincos_approximation(int neg, float128 r, uint64_t quotient)
192: {
193: if (quotient & 0x1) {
194: r = poly_cos(r);
195: neg = 0;
196: } else {
197: r = poly_sin(r);
198: }
199:
200: floatx80 result = float128_to_floatx80(r);
201: if (quotient & 0x2)
202: neg = ! neg;
203:
204: if (neg)
205: result = floatx80_chs(result);
206:
207: return result;
208: }
209:
210: // =================================================
211: // SFFSINCOS Compute sin(x) and cos(x)
212: // =================================================
213:
214: //
215: // Uses the following identities:
216: // ----------------------------------------------------------
217: //
218: // sin(-x) = -sin(x)
219: // cos(-x) = cos(x)
220: //
221: // sin(x+y) = sin(x)*cos(y)+cos(x)*sin(y)
222: // cos(x+y) = sin(x)*sin(y)+cos(x)*cos(y)
223: //
224: // sin(x+ pi/2) = cos(x)
225: // sin(x+ pi) = -sin(x)
226: // sin(x+3pi/2) = -cos(x)
227: // sin(x+2pi) = sin(x)
228: //
229:
230: int sf_fsincos(floatx80 a, floatx80 *sin_a, floatx80 *cos_a)
231: {
232: uint64_t aSig0, aSig1 = 0;
233: int32_t aExp, zExp, expDiff;
234: int aSign, zSign;
235: int q = 0;
236:
237: aSig0 = extractFloatx80Frac(a);
238: aExp = extractFloatx80Exp(a);
239: aSign = extractFloatx80Sign(a);
240:
241: /* invalid argument */
242: if (aExp == 0x7FFF) {
243: if ((uint64_t) (aSig0<<1)) {
244: sincos_invalid(sin_a, cos_a, propagateFloatx80NaNOneArg(a));
245: return 0;
246: }
247:
248: float_raise(float_flag_invalid);
249: sincos_invalid(sin_a, cos_a, floatx80_default_nan);
250: return 0;
251: }
252:
253: if (aExp == 0) {
254: if (aSig0 == 0) {
255: sincos_tiny_argument(sin_a, cos_a, a);
256: return 0;
257: }
258:
259: // float_raise(float_flag_denormal);
260:
261: /* handle pseudo denormals */
262: if (! (aSig0 & 0x8000000000000000U))
263: {
264: float_raise(float_flag_inexact);
265: if (sin_a)
266: float_raise(float_flag_underflow);
267: sincos_tiny_argument(sin_a, cos_a, a);
268: return 0;
269: }
270:
271: normalizeFloatx80Subnormal(aSig0, &aExp, &aSig0);
272: }
273:
274: zSign = aSign;
275: zExp = EXP_BIAS;
276: expDiff = aExp - zExp;
277:
278: /* argument is out-of-range */
279: if (expDiff >= 63)
280: return -1;
281:
282: float_raise(float_flag_inexact);
283:
284: if (expDiff < -1) { // doesn't require reduction
285: if (expDiff <= -68) {
286: a = packFloatx80(aSign, aExp, aSig0);
287: sincos_tiny_argument(sin_a, cos_a, a);
288: return 0;
289: }
290: zExp = aExp;
291: }
292: else {
293: q = reduce_trig_arg(expDiff, &zSign, &aSig0, &aSig1);
294: }
295:
296: /* **************************** */
297: /* argument reduction completed */
298: /* **************************** */
299:
300: /* using float128 for approximation */
301: float128 r = normalizeRoundAndPackFloat128(0, zExp-0x10, aSig0, aSig1);
302:
303: if (aSign) q = -q;
304: if (cos_a) *cos_a = sincos_approximation(zSign, r, q+1);
305: if (sin_a) *sin_a = sincos_approximation(zSign, r, q);
306:
307: return 0;
308: }
309:
310: int floatx80_fsincos(floatx80 a, floatx80 *sin_a, floatx80 *cos_a)
311: {
312: return sf_fsincos(a, sin_a, cos_a);
313: }
314:
315: int floatx80_fsin(floatx80 *a)
316: {
317: return sf_fsincos(*a, a, 0);
318: }
319:
320: int floatx80_fcos(floatx80 *a)
321: {
322: return sf_fsincos(*a, 0, a);
323: }
324:
325: // =================================================
326: // FPTAN Compute tan(x)
327: // =================================================
328:
329: //
330: // Uses the following identities:
331: //
332: // 1. ----------------------------------------------------------
333: //
334: // sin(-x) = -sin(x)
335: // cos(-x) = cos(x)
336: //
337: // sin(x+y) = sin(x)*cos(y)+cos(x)*sin(y)
338: // cos(x+y) = sin(x)*sin(y)+cos(x)*cos(y)
339: //
340: // sin(x+ pi/2) = cos(x)
341: // sin(x+ pi) = -sin(x)
342: // sin(x+3pi/2) = -cos(x)
343: // sin(x+2pi) = sin(x)
344: //
345: // 2. ----------------------------------------------------------
346: //
347: // sin(x)
348: // tan(x) = ------
349: // cos(x)
350: //
351:
352: int floatx80_ftan(floatx80 *a)
353: {
354: uint64_t aSig0, aSig1 = 0;
355: int32_t aExp, zExp, expDiff;
356: int aSign, zSign;
357: int q = 0;
358:
359: aSig0 = extractFloatx80Frac(*a);
360: aExp = extractFloatx80Exp(*a);
361: aSign = extractFloatx80Sign(*a);
362:
363: /* invalid argument */
364: if (aExp == 0x7FFF) {
365: if ((uint64_t) (aSig0<<1))
366: {
367: *a = propagateFloatx80NaNOneArg(*a);
368: return 0;
369: }
370:
371: float_raise(float_flag_invalid);
372: *a = floatx80_default_nan;
373: return 0;
374: }
375:
376: if (aExp == 0) {
377: if (aSig0 == 0) return 0;
378: // float_raise(float_flag_denormal);
379: /* handle pseudo denormals */
380: if (! (aSig0 & 0x8000000000000000U))
381: {
382: float_raise(float_flag_inexact | float_flag_underflow);
383: return 0;
384: }
385: normalizeFloatx80Subnormal(aSig0, &aExp, &aSig0);
386: }
387:
388: zSign = aSign;
389: zExp = EXP_BIAS;
390: expDiff = aExp - zExp;
391:
392: /* argument is out-of-range */
393: if (expDiff >= 63)
394: return -1;
395:
396: float_raise(float_flag_inexact);
397:
398: if (expDiff < -1) { // doesn't require reduction
399: if (expDiff <= -68) {
400: *a = packFloatx80(aSign, aExp, aSig0);
401: return 0;
402: }
403: zExp = aExp;
404: }
405: else {
406: q = reduce_trig_arg(expDiff, &zSign, &aSig0, &aSig1);
407: }
408:
409: /* **************************** */
410: /* argument reduction completed */
411: /* **************************** */
412:
413: /* using float128 for approximation */
414: float128 r = normalizeRoundAndPackFloat128(0, zExp-0x10, aSig0, aSig1);
415:
416: float128 sin_r = poly_sin(r);
417: float128 cos_r = poly_cos(r);
418:
419: if (q & 0x1) {
420: r = float128_div(cos_r, sin_r);
421: zSign = ! zSign;
422: } else {
423: r = float128_div(sin_r, cos_r);
424: }
425:
426: *a = float128_to_floatx80(r);
427: if (zSign)
428: *a = floatx80_chs(*a);
429:
430: return 0;
431: }
432:
433: // 2 3 4 n
434: // f(x) ~ C + (C * x) + (C * x) + (C * x) + (C * x) + ... + (C * x)
435: // 0 1 2 3 4 n
436: //
437: // -- 2k -- 2k+1
438: // p(x) = > C * x q(x) = > C * x
439: // -- 2k -- 2k+1
440: //
441: // f(x) ~ [ p(x) + x * q(x) ]
442: //
443:
444: float128 EvalPoly(float128 x, float128 *arr, unsigned n)
445: {
446: float128 x2 = float128_mul(x, x);
447: unsigned i;
448:
449: assert(n > 1);
450:
451: float128 r1 = arr[--n];
452: i = n;
453: while(i >= 2) {
454: r1 = float128_mul(r1, x2);
455: i -= 2;
456: r1 = float128_add(r1, arr[i]);
457: }
458: if (i) r1 = float128_mul(r1, x);
459:
460: float128 r2 = arr[--n];
461: i = n;
462: while(i >= 2) {
463: r2 = float128_mul(r2, x2);
464: i -= 2;
465: r2 = float128_add(r2, arr[i]);
466: }
467: if (i) r2 = float128_mul(r2, x);
468:
469: return float128_add(r1, r2);
470: }
471:
472: // 2 4 6 8 2n
473: // f(x) ~ C + (C * x) + (C * x) + (C * x) + (C * x) + ... + (C * x)
474: // 0 1 2 3 4 n
475: //
476: // -- 4k -- 4k+2
477: // p(x) = > C * x q(x) = > C * x
478: // -- 2k -- 2k+1
479: //
480: // 2
481: // f(x) ~ [ p(x) + x * q(x) ]
482: //
483:
484: float128 EvenPoly(float128 x, float128 *arr, unsigned n)
485: {
486: return EvalPoly(float128_mul(x, x), arr, n);
487: }
488:
489: // 3 5 7 9 2n+1
490: // f(x) ~ (C * x) + (C * x) + (C * x) + (C * x) + (C * x) + ... + (C * x)
491: // 0 1 2 3 4 n
492: // 2 4 6 8 2n
493: // = x * [ C + (C * x) + (C * x) + (C * x) + (C * x) + ... + (C * x)
494: // 0 1 2 3 4 n
495: //
496: // -- 4k -- 4k+2
497: // p(x) = > C * x q(x) = > C * x
498: // -- 2k -- 2k+1
499: //
500: // 2
501: // f(x) ~ x * [ p(x) + x * q(x) ]
502: //
503:
504: float128 OddPoly(float128 x, float128 *arr, unsigned n)
505: {
506: return float128_mul(x, EvenPoly(x, arr, n));
507: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.