![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to softfloat_fpsp.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 softfloat_fpsp.c CVS log | Up to [Previous NeXT emulator] / previous / src / softfloat |
1.1 root 1:
2: /*============================================================================
3:
4: This C source file is an extension to the SoftFloat IEC/IEEE Floating-point
5: Arithmetic Package, Release 2a.
6:
7: Written by Andreas Grabher for Previous, NeXT Computer Emulator.
8:
9: =============================================================================*/
10:
11: #include "softfloat.h"
12: #include "softfloat_fpsp_tables.h"
13:
14:
15: /*----------------------------------------------------------------------------
16: | Algorithms for transcendental functions supported by MC68881 and MC68882
17: | mathematical coprocessors. The functions are derived from FPSP library.
18: *----------------------------------------------------------------------------*/
19:
20: #define pi_sig LIT64(0xc90fdaa22168c235)
21: #define pi_sig0 LIT64(0xc90fdaa22168c234)
22: #define pi_sig1 LIT64(0xc4c6628b80dc1cd1)
23:
24: #define pi_exp 0x4000
25: #define piby2_exp 0x3FFF
26: #define piby4_exp 0x3FFE
27:
28: #define one_exp 0x3FFF
29: #define one_sig LIT64(0x8000000000000000)
30:
31:
32: /*----------------------------------------------------------------------------
33: | Function for compactifying extended double-precision floating point values.
34: *----------------------------------------------------------------------------*/
35:
36: int32 floatx80_make_compact(int32 aExp, bits64 aSig)
37: {
38: return (aExp<<16)|(aSig>>48);
39: }
40:
41:
42: /*----------------------------------------------------------------------------
43: | Arc cosine
44: *----------------------------------------------------------------------------*/
45:
46: floatx80 floatx80_acos(floatx80 a)
47: {
48: flag aSign;
49: int32 aExp;
50: bits64 aSig;
51:
52: int8 user_rnd_mode, user_rnd_prec;
53:
54: int32 compact;
55: floatx80 fp0, fp1, one;
56:
57: aSig = extractFloatx80Frac(a);
58: aExp = extractFloatx80Exp(a);
59: aSign = extractFloatx80Sign(a);
60:
61: if (aExp == 0x7FFF && (bits64) (aSig<<1)) {
62: return propagateFloatx80NaNOneArg(a);
63: }
64: if (aExp == 0 && aSig == 0) {
65: float_raise(float_flag_inexact);
66: return roundAndPackFloatx80(floatx80_rounding_precision, 0, piby2_exp, pi_sig, 0);
67: }
68:
69: compact = floatx80_make_compact(aExp, aSig);
70:
71: if (compact >= 0x3FFF8000) { // |X| >= 1
72: if (aExp == one_exp && aSig == one_sig) { // |X| == 1
73: if (aSign) { // X == -1
74: a = packFloatx80(0, pi_exp, pi_sig);
75: float_raise(float_flag_inexact);
76: return floatx80_move(a);
77: } else { // X == +1
78: return packFloatx80(0, 0, 0);
79: }
80: } else { // |X| > 1
81: float_raise(float_flag_invalid);
82: a.low = floatx80_default_nan_low;
83: a.high = floatx80_default_nan_high;
84: return a;
85: }
86: } // |X| < 1
87:
88: user_rnd_mode = float_rounding_mode;
89: user_rnd_prec = floatx80_rounding_precision;
90: float_rounding_mode = float_round_nearest_even;
91: floatx80_rounding_precision = 80;
92:
93: one = packFloatx80(0, one_exp, one_sig);
94: fp0 = a;
95:
96: fp1 = floatx80_add(one, fp0); // 1 + X
97: fp0 = floatx80_sub(one, fp0); // 1 - X
98: fp0 = floatx80_div(fp0, fp1); // (1-X)/(1+X)
99: fp0 = floatx80_sqrt(fp0); // SQRT((1-X)/(1+X))
100: fp0 = floatx80_atan(fp0); // ATAN(SQRT((1-X)/(1+X)))
101:
102: float_rounding_mode = user_rnd_mode;
103: floatx80_rounding_precision = user_rnd_prec;
104:
105: a = floatx80_add(fp0, fp0); // 2 * ATAN(SQRT((1-X)/(1+X)))
106:
107: float_raise(float_flag_inexact);
108:
109: return a;
110: }
111:
112: /*----------------------------------------------------------------------------
113: | Arc sine
114: *----------------------------------------------------------------------------*/
115:
116: floatx80 floatx80_asin(floatx80 a)
117: {
118: flag aSign;
119: int32 aExp;
120: bits64 aSig;
121:
122: int8 user_rnd_mode, user_rnd_prec;
123:
124: int32 compact;
125: floatx80 fp0, fp1, fp2, one;
126:
127: aSig = extractFloatx80Frac(a);
128: aExp = extractFloatx80Exp(a);
129: aSign = extractFloatx80Sign(a);
130:
131: if (aExp == 0x7FFF && (bits64) (aSig<<1)) {
132: return propagateFloatx80NaNOneArg(a);
133: }
134:
135: if (aExp == 0 && aSig == 0) {
136: return packFloatx80(aSign, 0, 0);
137: }
138:
139: compact = floatx80_make_compact(aExp, aSig);
140:
141: if (compact >= 0x3FFF8000) { // |X| >= 1
142: if (aExp == one_exp && aSig == one_sig) { // |X| == 1
143: float_raise(float_flag_inexact);
144: a = packFloatx80(aSign, piby2_exp, pi_sig);
145: return floatx80_move(a);
146: } else { // |X| > 1
147: float_raise(float_flag_invalid);
148: a.low = floatx80_default_nan_low;
149: a.high = floatx80_default_nan_high;
150: return a;
151: }
152:
153: } // |X| < 1
154:
155: user_rnd_mode = float_rounding_mode;
156: user_rnd_prec = floatx80_rounding_precision;
157: float_rounding_mode = float_round_nearest_even;
158: floatx80_rounding_precision = 80;
159:
160: one = packFloatx80(0, one_exp, one_sig);
161: fp0 = a;
162:
163: fp1 = floatx80_sub(one, fp0); // 1 - X
164: fp2 = floatx80_add(one, fp0); // 1 + X
165: fp1 = floatx80_mul(fp2, fp1); // (1+X)*(1-X)
166: fp1 = floatx80_sqrt(fp1); // SQRT((1+X)*(1-X))
167: fp0 = floatx80_div(fp0, fp1); // X/SQRT((1+X)*(1-X))
168:
169: float_rounding_mode = user_rnd_mode;
170: floatx80_rounding_precision = user_rnd_prec;
171:
172: a = floatx80_atan(fp0); // ATAN(X/SQRT((1+X)*(1-X)))
173:
174: float_raise(float_flag_inexact);
175:
176: return a;
177: }
178:
179: /*----------------------------------------------------------------------------
180: | Arc tangent
181: *----------------------------------------------------------------------------*/
182:
183: floatx80 floatx80_atan(floatx80 a)
184: {
185: flag aSign;
186: int32 aExp;
187: bits64 aSig;
188:
189: int8 user_rnd_mode, user_rnd_prec;
190:
191: int32 compact, tbl_index;
192: floatx80 fp0, fp1, fp2, fp3, xsave;
193:
194: aSig = extractFloatx80Frac(a);
195: aExp = extractFloatx80Exp(a);
196: aSign = extractFloatx80Sign(a);
197:
198: if (aExp == 0x7FFF) {
199: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a);
200: a = packFloatx80(aSign, piby2_exp, pi_sig);
201: float_raise(float_flag_inexact);
202: return floatx80_move(a);
203: }
204:
205: if (aExp == 0 && aSig == 0) {
206: return packFloatx80(aSign, 0, 0);
207: }
208:
209: compact = floatx80_make_compact(aExp, aSig);
210:
211: user_rnd_mode = float_rounding_mode;
212: user_rnd_prec = floatx80_rounding_precision;
213: float_rounding_mode = float_round_nearest_even;
214: floatx80_rounding_precision = 80;
215:
216: if (compact < 0x3FFB8000 || compact > 0x4002FFFF) { // |X| >= 16 or |X| < 1/16
217: if (compact > 0x3FFF8000) { // |X| >= 16
218: if (compact > 0x40638000) { // |X| > 2^(100)
219: fp0 = packFloatx80(aSign, piby2_exp, pi_sig);
220: fp1 = packFloatx80(aSign, 0x0001, one_sig);
221:
222: float_rounding_mode = user_rnd_mode;
223: floatx80_rounding_precision = user_rnd_prec;
224:
225: a = floatx80_sub(fp0, fp1);
226:
227: float_raise(float_flag_inexact);
228:
229: return a;
230: } else {
231: fp0 = a;
232: fp1 = packFloatx80(1, one_exp, one_sig); // -1
233: fp1 = floatx80_div(fp1, fp0); // X' = -1/X
234: xsave = fp1;
235: fp0 = floatx80_mul(fp1, fp1); // Y = X'*X'
236: fp1 = floatx80_mul(fp0, fp0); // Z = Y*Y
237: fp3 = float64_to_floatx80(LIT64(0xBFB70BF398539E6A)); // C5
238: fp2 = float64_to_floatx80(LIT64(0x3FBC7187962D1D7D)); // C4
239: fp3 = floatx80_mul(fp3, fp1); // Z*C5
240: fp2 = floatx80_mul(fp2, fp1); // Z*C4
241: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBFC24924827107B8))); // C3+Z*C5
242: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FC999999996263E))); // C2+Z*C4
243: fp1 = floatx80_mul(fp1, fp3); // Z*(C3+Z*C5)
244: fp2 = floatx80_mul(fp2, fp0); // Y*(C2+Z*C4)
245: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0xBFD5555555555536))); // C1+Z*(C3+Z*C5)
246: fp0 = floatx80_mul(fp0, xsave); // X'*Y
247: fp1 = floatx80_add(fp1, fp2); // [Y*(C2+Z*C4)]+[C1+Z*(C3+Z*C5)]
248: fp0 = floatx80_mul(fp0, fp1); // X'*Y*([B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))]) ??
249: fp0 = floatx80_add(fp0, xsave);
250: fp1 = packFloatx80(aSign, piby2_exp, pi_sig);
251:
252: float_rounding_mode = user_rnd_mode;
253: floatx80_rounding_precision = user_rnd_prec;
254:
255: a = floatx80_add(fp0, fp1);
256:
257: float_raise(float_flag_inexact);
258:
259: return a;
260: }
261: } else { // |X| < 1/16
262: if (compact < 0x3FD78000) { // |X| < 2^(-40)
263: float_rounding_mode = user_rnd_mode;
264: floatx80_rounding_precision = user_rnd_prec;
265:
266: a = floatx80_move(a);
267:
268: float_raise(float_flag_inexact);
269:
270: return a;
271: } else {
272: fp0 = a;
273: xsave = a;
274: fp0 = floatx80_mul(fp0, fp0); // Y = X*X
275: fp1 = floatx80_mul(fp0, fp0); // Z = Y*Y
276: fp2 = float64_to_floatx80(LIT64(0x3FB344447F876989)); // B6
277: fp3 = float64_to_floatx80(LIT64(0xBFB744EE7FAF45DB)); // B5
278: fp2 = floatx80_mul(fp2, fp1); // Z*B6
279: fp3 = floatx80_mul(fp3, fp1); // Z*B5
280: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FBC71C646940220))); // B4+Z*B6
281: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBFC24924921872F9))); // B3+Z*B5
282: fp2 = floatx80_mul(fp2, fp1); // Z*(B4+Z*B6)
283: fp1 = floatx80_mul(fp1, fp3); // Z*(B3+Z*B5)
284: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FC9999999998FA9))); // B2+Z*(B4+Z*B6)
285: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0xBFD5555555555555))); // B1+Z*(B3+Z*B5)
286: fp2 = floatx80_mul(fp2, fp0); // Y*(B2+Z*(B4+Z*B6))
287: fp0 = floatx80_mul(fp0, xsave); // X*Y
288: fp1 = floatx80_add(fp1, fp2); // [B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))]
289: fp0 = floatx80_mul(fp0, fp1); // X*Y*([B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))])
290:
291: float_rounding_mode = user_rnd_mode;
292: floatx80_rounding_precision = user_rnd_prec;
293:
294: a = floatx80_add(fp0, xsave);
295:
296: float_raise(float_flag_inexact);
297:
298: return a;
299: }
300: }
301: } else {
302: aSig &= LIT64(0xF800000000000000);
303: aSig |= LIT64(0x0400000000000000);
304: xsave = packFloatx80(aSign, aExp, aSig); // F
305: fp0 = a;
306: fp1 = a; // X
307: fp2 = packFloatx80(0, one_exp, one_sig); // 1
308: fp1 = floatx80_mul(fp1, xsave); // X*F
309: fp0 = floatx80_sub(fp0, xsave); // X-F
310: fp1 = floatx80_add(fp1, fp2); // 1 + X*F
311: fp0 = floatx80_div(fp0, fp1); // U = (X-F)/(1+X*F)
312:
313: tbl_index = compact;
314:
315: tbl_index &= 0x7FFF0000;
316: tbl_index -= 0x3FFB0000;
317: tbl_index >>= 1;
318: tbl_index += compact&0x00007800;
319: tbl_index >>= 11;
320:
321: fp3 = atan_tbl[tbl_index];
322:
323: fp3.high |= aSign ? 0x8000 : 0; // ATAN(F)
324:
325: fp1 = floatx80_mul(fp0, fp0); // V = U*U
326: fp2 = float64_to_floatx80(LIT64(0xBFF6687E314987D8)); // A3
327: fp2 = floatx80_add(fp2, fp1); // A3+V
328: fp2 = floatx80_mul(fp2, fp1); // V*(A3+V)
329: fp1 = floatx80_mul(fp1, fp0); // U*V
330: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x4002AC6934A26DB3))); // A2+V*(A3+V)
331: fp1 = floatx80_mul(fp1, float64_to_floatx80(LIT64(0xBFC2476F4E1DA28E))); // A1+U*V
332: fp1 = floatx80_mul(fp1, fp2); // A1*U*V*(A2+V*(A3+V))
333: fp0 = floatx80_add(fp0, fp1); // ATAN(U)
334:
335: float_rounding_mode = user_rnd_mode;
336: floatx80_rounding_precision = user_rnd_prec;
337:
338: a = floatx80_add(fp0, fp3); // ATAN(X)
339:
340: float_raise(float_flag_inexact);
341:
342: return a;
343: }
344: }
345:
346: /*----------------------------------------------------------------------------
347: | Hyperbolic arc tangent
348: *----------------------------------------------------------------------------*/
349:
350: floatx80 floatx80_atanh(floatx80 a)
351: {
352: flag aSign;
353: int32 aExp;
354: bits64 aSig;
355:
356: int8 user_rnd_mode, user_rnd_prec;
357:
358: int32 compact;
359: floatx80 fp0, fp1, fp2, one;
360:
361: aSig = extractFloatx80Frac(a);
362: aExp = extractFloatx80Exp(a);
363: aSign = extractFloatx80Sign(a);
364:
365: if (aExp == 0x7FFF && (bits64) (aSig<<1)) {
366: return propagateFloatx80NaNOneArg(a);
367: }
368:
369: if (aExp == 0 && aSig == 0) {
370: return packFloatx80(aSign, 0, 0);
371: }
372:
373: compact = floatx80_make_compact(aExp, aSig);
374:
375: if (compact >= 0x3FFF8000) { // |X| >= 1
376: if (aExp == one_exp && aSig == one_sig) { // |X| == 1
377: float_raise(float_flag_divbyzero);
378: return packFloatx80(aSign, 0x7FFF, floatx80_default_infinity_low);
379: } else { // |X| > 1
380: float_raise(float_flag_invalid);
381: a.low = floatx80_default_nan_low;
382: a.high = floatx80_default_nan_high;
383: return a;
384: }
385: } // |X| < 1
386:
387: user_rnd_mode = float_rounding_mode;
388: user_rnd_prec = floatx80_rounding_precision;
389: float_rounding_mode = float_round_nearest_even;
390: floatx80_rounding_precision = 80;
391:
392: one = packFloatx80(0, one_exp, one_sig);
393: fp2 = packFloatx80(aSign, 0x3FFE, one_sig); // SIGN(X) * (1/2)
394: fp0 = packFloatx80(0, aExp, aSig); // Y = |X|
395: fp1 = packFloatx80(1, aExp, aSig); // -Y
396: fp0 = floatx80_add(fp0, fp0); // 2Y
397: fp1 = floatx80_add(fp1, one); // 1-Y
398: fp0 = floatx80_div(fp0, fp1); // Z = 2Y/(1-Y)
399: fp0 = floatx80_lognp1(fp0); // LOG1P(Z)
400:
401: float_rounding_mode = user_rnd_mode;
402: floatx80_rounding_precision = user_rnd_prec;
403:
404: a = floatx80_mul(fp0, fp2); // ATANH(X) = SIGN(X) * (1/2) * LOG1P(Z)
405:
406: float_raise(float_flag_inexact);
407:
408: return a;
409: }
410:
411: /*----------------------------------------------------------------------------
412: | Cosine
413: *----------------------------------------------------------------------------*/
414:
415: floatx80 floatx80_cos(floatx80 a)
416: {
417: flag aSign, xSign;
418: int32 aExp, xExp;
419: bits64 aSig, xSig;
420:
421: int8 user_rnd_mode, user_rnd_prec;
422:
423: int32 compact, l, n, j;
424: floatx80 fp0, fp1, fp2, fp3, fp4, fp5, x, invtwopi, twopi1, twopi2;
425: float32 posneg1, twoto63;
426: flag adjn, endflag;
427:
428: aSig = extractFloatx80Frac(a);
429: aExp = extractFloatx80Exp(a);
430: aSign = extractFloatx80Sign(a);
431:
432: if (aExp == 0x7FFF) {
433: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a);
434: float_raise(float_flag_invalid);
435: a.low = floatx80_default_nan_low;
436: a.high = floatx80_default_nan_high;
437: return a;
438: }
439:
440: if (aExp == 0 && aSig == 0) {
441: return packFloatx80(0, one_exp, one_sig);
442: }
443:
444: adjn = 1;
445:
446: user_rnd_mode = float_rounding_mode;
447: user_rnd_prec = floatx80_rounding_precision;
448: float_rounding_mode = float_round_nearest_even;
449: floatx80_rounding_precision = 80;
450:
451: compact = floatx80_make_compact(aExp, aSig);
452:
453: fp0 = a;
454:
455: if (compact < 0x3FD78000 || compact > 0x4004BC7E) { // 2^(-40) > |X| > 15 PI
456: if (compact > 0x3FFF8000) { // |X| >= 15 PI
457: // REDUCEX
458: fp1 = packFloatx80(0, 0, 0);
459: if (compact == 0x7FFEFFFF) {
460: twopi1 = packFloatx80(aSign ^ 1, 0x7FFE, LIT64(0xC90FDAA200000000));
461: twopi2 = packFloatx80(aSign ^ 1, 0x7FDC, LIT64(0x85A308D300000000));
462: fp0 = floatx80_add(fp0, twopi1);
463: fp1 = fp0;
464: fp0 = floatx80_add(fp0, twopi2);
465: fp1 = floatx80_sub(fp1, fp0);
466: fp1 = floatx80_add(fp1, twopi2);
467: }
468: loop:
469: xSign = extractFloatx80Sign(fp0);
470: xExp = extractFloatx80Exp(fp0);
471: xExp -= 0x3FFF;
472: if (xExp <= 28) {
473: l = 0;
474: endflag = 1;
475: } else {
476: l = xExp - 27;
477: endflag = 0;
478: }
479: invtwopi = packFloatx80(0, 0x3FFE - l, LIT64(0xA2F9836E4E44152A)); // INVTWOPI
480: twopi1 = packFloatx80(0, 0x3FFF + l, LIT64(0xC90FDAA200000000));
481: twopi2 = packFloatx80(0, 0x3FDD + l, LIT64(0x85A308D300000000));
482:
483: twoto63 = 0x5F000000;
484: twoto63 |= xSign ? 0x80000000 : 0x00000000; // SIGN(INARG)*2^63 IN SGL
485:
486: fp2 = floatx80_mul(fp0, invtwopi);
487: fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63)); // THE FRACTIONAL PART OF FP2 IS ROUNDED
488: fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63)); // FP2 is N
489: fp4 = floatx80_mul(twopi1, fp2); // W = N*P1
490: fp5 = floatx80_mul(twopi2, fp2); // w = N*P2
491: fp3 = floatx80_add(fp4, fp5); // FP3 is P
492: fp4 = floatx80_sub(fp4, fp3); // W-P
493: fp0 = floatx80_sub(fp0, fp3); // FP0 is A := R - P
494: fp4 = floatx80_add(fp4, fp5); // FP4 is p = (W-P)+w
495: fp3 = fp0; // FP3 is A
496: fp1 = floatx80_sub(fp1, fp4); // FP1 is a := r - p
497: fp0 = floatx80_add(fp0, fp1); // FP0 is R := A+a
498:
499: if (endflag > 0) {
500: n = floatx80_to_int32(fp2);
501: goto sincont;
502: }
503: fp3 = floatx80_sub(fp3, fp0); // A-R
504: fp1 = floatx80_add(fp1, fp3); // FP1 is r := (A-R)+a
505: goto loop;
506: } else {
507: // SINSM
508: fp0 = float32_to_floatx80(0x3F800000); // 1
509:
510: float_rounding_mode = user_rnd_mode;
511: floatx80_rounding_precision = user_rnd_prec;
512:
513: if (adjn) {
514: // COSTINY
515: a = floatx80_sub(fp0, float32_to_floatx80(0x00800000));
516: } else {
517: // SINTINY
518: a = floatx80_move(a);
519: }
520: float_raise(float_flag_inexact);
521:
522: return a;
523: }
524: } else {
525: fp1 = floatx80_mul(fp0, float64_to_floatx80(LIT64(0x3FE45F306DC9C883))); // X*2/PI
526:
527: n = floatx80_to_int32(fp1);
528: j = 32 + n;
529:
530: fp0 = floatx80_sub(fp0, pi_tbl[j]); // X-Y1
531: fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j])); // FP0 IS R = (X-Y1)-Y2
532:
533: sincont:
534: if ((n + adjn) & 1) {
535: // COSPOLY
536: fp0 = floatx80_mul(fp0, fp0); // FP0 IS S
537: fp1 = floatx80_mul(fp0, fp0); // FP1 IS T
538: fp2 = float64_to_floatx80(LIT64(0x3D2AC4D0D6011EE3)); // B8
539: fp3 = float64_to_floatx80(LIT64(0xBDA9396F9F45AC19)); // B7
540:
541: xSign = extractFloatx80Sign(fp0); // X IS S
542: xExp = extractFloatx80Exp(fp0);
543: xSig = extractFloatx80Frac(fp0);
544:
545: if (((n + adjn) >> 1) & 1) {
546: xSign ^= 1;
547: posneg1 = 0xBF800000; // -1
548: } else {
549: xSign ^= 0;
550: posneg1 = 0x3F800000; // 1
551: } // X IS NOW R'= SGN*R
552:
553: fp2 = floatx80_mul(fp2, fp1); // TB8
554: fp3 = floatx80_mul(fp3, fp1); // TB7
555: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3E21EED90612C972))); // B6+TB8
556: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBE927E4FB79D9FCF))); // B5+TB7
557: fp2 = floatx80_mul(fp2, fp1); // T(B6+TB8)
558: fp3 = floatx80_mul(fp3, fp1); // T(B5+TB7)
559: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3EFA01A01A01D423))); // B4+T(B6+TB8)
560: fp4 = packFloatx80(1, 0x3FF5, LIT64(0xB60B60B60B61D438));
561: fp3 = floatx80_add(fp3, fp4); // B3+T(B5+TB7)
562: fp2 = floatx80_mul(fp2, fp1); // T(B4+T(B6+TB8))
563: fp1 = floatx80_mul(fp1, fp3); // T(B3+T(B5+TB7))
564: fp4 = packFloatx80(0, 0x3FFA, LIT64(0xAAAAAAAAAAAAAB5E));
565: fp2 = floatx80_add(fp2, fp4); // B2+T(B4+T(B6+TB8))
566: fp1 = floatx80_add(fp1, float32_to_floatx80(0xBF000000)); // B1+T(B3+T(B5+TB7))
567: fp0 = floatx80_mul(fp0, fp2); // S(B2+T(B4+T(B6+TB8)))
568: fp0 = floatx80_add(fp0, fp1); // [B1+T(B3+T(B5+TB7))]+[S(B2+T(B4+T(B6+TB8)))]
569:
570: x = packFloatx80(xSign, xExp, xSig);
571: fp0 = floatx80_mul(fp0, x);
572:
573: float_rounding_mode = user_rnd_mode;
574: floatx80_rounding_precision = user_rnd_prec;
575:
576: a = floatx80_add(fp0, float32_to_floatx80(posneg1));
577:
578: float_raise(float_flag_inexact);
579:
580: return a;
581: } else {
582: // SINPOLY
583: xSign = extractFloatx80Sign(fp0); // X IS R
584: xExp = extractFloatx80Exp(fp0);
585: xSig = extractFloatx80Frac(fp0);
586:
587: xSign ^= ((n + adjn) >> 1) & 1; // X IS NOW R'= SGN*R
588:
589: fp0 = floatx80_mul(fp0, fp0); // FP0 IS S
590: fp1 = floatx80_mul(fp0, fp0); // FP1 IS T
591: fp3 = float64_to_floatx80(LIT64(0xBD6AAA77CCC994F5)); // A7
592: fp2 = float64_to_floatx80(LIT64(0x3DE612097AAE8DA1)); // A6
593: fp3 = floatx80_mul(fp3, fp1); // T*A7
594: fp2 = floatx80_mul(fp2, fp1); // T*A6
595: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBE5AE6452A118AE4))); // A5+T*A7
596: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3EC71DE3A5341531))); // A4+T*A6
597: fp3 = floatx80_mul(fp3, fp1); // T(A5+TA7)
598: fp2 = floatx80_mul(fp2, fp1); // T(A4+TA6)
599: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBF2A01A01A018B59))); // A3+T(A5+TA7)
600: fp4 = packFloatx80(0, 0x3FF8, LIT64(0x88888888888859AF));
601: fp2 = floatx80_add(fp2, fp4); // A2+T(A4+TA6)
602: fp1 = floatx80_mul(fp1, fp3); // T(A3+T(A5+TA7))
603: fp2 = floatx80_mul(fp2, fp0); // S(A2+T(A4+TA6))
604: fp4 = packFloatx80(1, 0x3FFC, LIT64(0xAAAAAAAAAAAAAA99));
605: fp1 = floatx80_add(fp1, fp4); // A1+T(A3+T(A5+TA7))
606: fp1 = floatx80_add(fp1, fp2); // [A1+T(A3+T(A5+TA7))]+[S(A2+T(A4+TA6))]
607:
608: x = packFloatx80(xSign, xExp, xSig);
609: fp0 = floatx80_mul(fp0, x); // R'*S
610: fp0 = floatx80_mul(fp0, fp1); // SIN(R')-R'
611:
612: float_rounding_mode = user_rnd_mode;
613: floatx80_rounding_precision = user_rnd_prec;
614:
615: a = floatx80_add(fp0, x);
616:
617: float_raise(float_flag_inexact);
618:
619: return a;
620: }
621: }
622: }
623:
624: /*----------------------------------------------------------------------------
625: | Hyperbolic cosine
626: *----------------------------------------------------------------------------*/
627:
628: floatx80 floatx80_cosh(floatx80 a)
629: {
630: flag aSign;
631: int32 aExp;
632: bits64 aSig;
633:
634: int8 user_rnd_mode, user_rnd_prec;
635:
636: int32 compact;
637: floatx80 fp0, fp1;
638:
639: aSig = extractFloatx80Frac(a);
640: aExp = extractFloatx80Exp(a);
641: aSign = extractFloatx80Sign(a);
642:
643: if (aExp == 0x7FFF) {
644: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a);
645: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low);
646: }
647:
648: if (aExp == 0 && aSig == 0) {
649: return packFloatx80(0, one_exp, one_sig);
650: }
651:
652: user_rnd_mode = float_rounding_mode;
653: user_rnd_prec = floatx80_rounding_precision;
654: float_rounding_mode = float_round_nearest_even;
655: floatx80_rounding_precision = 80;
656:
657: compact = floatx80_make_compact(aExp, aSig);
658:
659: if (compact > 0x400CB167) {
660: if (compact > 0x400CB2B3) {
661: float_rounding_mode = user_rnd_mode;
662: floatx80_rounding_precision = user_rnd_prec;
663: return roundAndPackFloatx80(floatx80_rounding_precision, 0, 0x8000, one_sig, 0);
664: } else {
665: fp0 = packFloatx80(0, aExp, aSig);
666: fp0 = floatx80_sub(fp0, float64_to_floatx80(LIT64(0x40C62D38D3D64634)));
667: fp0 = floatx80_sub(fp0, float64_to_floatx80(LIT64(0x3D6F90AEB1E75CC7)));
668: fp0 = floatx80_etox(fp0);
669: fp1 = packFloatx80(0, 0x7FFB, one_sig);
670:
671: float_rounding_mode = user_rnd_mode;
672: floatx80_rounding_precision = user_rnd_prec;
673:
674: a = floatx80_mul(fp0, fp1);
675:
676: float_raise(float_flag_inexact);
677:
678: return a;
679: }
680: }
681:
682: fp0 = packFloatx80(0, aExp, aSig); // |X|
683: fp0 = floatx80_etox(fp0); // EXP(|X|)
684: fp0 = floatx80_mul(fp0, float32_to_floatx80(0x3F000000)); // (1/2)*EXP(|X|)
685: fp1 = float32_to_floatx80(0x3E800000); // 1/4
686: fp1 = floatx80_div(fp1, fp0); // 1/(2*EXP(|X|))
687:
688: float_rounding_mode = user_rnd_mode;
689: floatx80_rounding_precision = user_rnd_prec;
690:
691: a = floatx80_add(fp0, fp1);
692:
693: float_raise(float_flag_inexact);
694:
695: return a;
696: }
697:
698: /*----------------------------------------------------------------------------
699: | e to x
700: *----------------------------------------------------------------------------*/
701:
702: floatx80 floatx80_etox(floatx80 a)
703: {
704: flag aSign;
705: int32 aExp;
706: bits64 aSig;
707:
708: int8 user_rnd_mode, user_rnd_prec;
709:
710: int32 compact, n, j, k, m, m1;
711: floatx80 fp0, fp1, fp2, fp3, l2, scale, adjscale;
712: flag adjflag;
713:
714: aSig = extractFloatx80Frac(a);
715: aExp = extractFloatx80Exp(a);
716: aSign = extractFloatx80Sign(a);
717:
718: if (aExp == 0x7FFF) {
719: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a);
720: if (aSign) return packFloatx80(0, 0, 0);
721: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low);
722: }
723:
724: if (aExp == 0 && aSig == 0) {
725: return packFloatx80(0, one_exp, one_sig);
726: }
727:
728: user_rnd_mode = float_rounding_mode;
729: user_rnd_prec = floatx80_rounding_precision;
730: float_rounding_mode = float_round_nearest_even;
731: floatx80_rounding_precision = 80;
732:
733: adjflag = 0;
734:
735: if (aExp >= 0x3FBE) { // |X| >= 2^(-65)
736: compact = floatx80_make_compact(aExp, aSig);
737:
738: if (compact < 0x400CB167) { // |X| < 16380 log2
739: fp0 = a;
740: fp1 = a;
741: fp0 = floatx80_mul(fp0, float32_to_floatx80(0x42B8AA3B)); // 64/log2 * X
742: adjflag = 0;
743: n = floatx80_to_int32(fp0); // int(64/log2*X)
744: fp0 = int32_to_floatx80(n);
745:
746: j = n & 0x3F; // J = N mod 64
747: m = n / 64; // NOTE: this is really arithmetic right shift by 6
748: if (n < 0 && j) { // arithmetic right shift is division and round towards minus infinity
749: m--;
750: }
751: m += 0x3FFF; // biased exponent of 2^(M)
752:
753: expcont1:
754: fp2 = fp0; // N
755: fp0 = floatx80_mul(fp0, float32_to_floatx80(0xBC317218)); // N * L1, L1 = lead(-log2/64)
756: l2 = packFloatx80(0, 0x3FDC, LIT64(0x82E308654361C4C6));
757: fp2 = floatx80_mul(fp2, l2); // N * L2, L1+L2 = -log2/64
758: fp0 = floatx80_add(fp0, fp1); // X + N*L1
759: fp0 = floatx80_add(fp0, fp2); // R
760:
761: fp1 = floatx80_mul(fp0, fp0); // S = R*R
762: fp2 = float32_to_floatx80(0x3AB60B70); // A5
763: fp2 = floatx80_mul(fp2, fp1); // fp2 is S*A5
764: fp3 = floatx80_mul(float32_to_floatx80(0x3C088895), fp1); // fp3 is S*A4
765: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FA5555555554431))); // fp2 is A3+S*A5
766: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0x3FC5555555554018))); // fp3 is A2+S*A4
767: fp2 = floatx80_mul(fp2, fp1); // fp2 is S*(A3+S*A5)
768: fp3 = floatx80_mul(fp3, fp1); // fp3 is S*(A2+S*A4)
769: fp2 = floatx80_add(fp2, float32_to_floatx80(0x3F000000)); // fp2 is A1+S*(A3+S*A5)
770: fp3 = floatx80_mul(fp3, fp0); // fp3 IS R*S*(A2+S*A4)
771: fp2 = floatx80_mul(fp2, fp1); // fp2 IS S*(A1+S*(A3+S*A5))
772: fp0 = floatx80_add(fp0, fp3); // fp0 IS R+R*S*(A2+S*A4)
773: fp0 = floatx80_add(fp0, fp2); // fp0 IS EXP(R) - 1
774:
775: fp1 = exp_tbl[j];
776: fp0 = floatx80_mul(fp0, fp1); // 2^(J/64)*(Exp(R)-1)
777: fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j])); // accurate 2^(J/64)
778: fp0 = floatx80_add(fp0, fp1); // 2^(J/64) + 2^(J/64)*(Exp(R)-1)
779:
780: scale = packFloatx80(0, m, one_sig);
781: if (adjflag) {
782: adjscale = packFloatx80(0, m1, one_sig);
783: fp0 = floatx80_mul(fp0, adjscale);
784: }
785:
786: float_rounding_mode = user_rnd_mode;
787: floatx80_rounding_precision = user_rnd_prec;
788:
789: a = floatx80_mul(fp0, scale);
790:
791: float_raise(float_flag_inexact);
792:
793: return a;
794: } else { // |X| >= 16380 log2
795: if (compact > 0x400CB27C) { // |X| >= 16480 log2
796: float_rounding_mode = user_rnd_mode;
797: floatx80_rounding_precision = user_rnd_prec;
798: if (aSign) {
799: a = roundAndPackFloatx80(floatx80_rounding_precision, 0, -0x1000, aSig, 0);
800: } else {
801: a = roundAndPackFloatx80(floatx80_rounding_precision, 0, 0x8000, aSig, 0);
802: }
803: float_raise(float_flag_inexact);
804:
805: return a;
806: } else {
807: fp0 = a;
808: fp1 = a;
809: fp0 = floatx80_mul(fp0, float32_to_floatx80(0x42B8AA3B)); // 64/log2 * X
810: adjflag = 1;
811: n = floatx80_to_int32(fp0); // int(64/log2*X)
812: fp0 = int32_to_floatx80(n);
813:
814: j = n & 0x3F; // J = N mod 64
815: k = n / 64; // NOTE: this is really arithmetic right shift by 6
816: if (n < 0 && j) { // arithmetic right shift is division and round towards minus infinity
817: k--;
818: }
819: m1 = k / 2; // NOTE: this is really arithmetic right shift by 1
820: if (k < 0 && (k & 1)) { // arithmetic right shift is division and round towards minus infinity
821: m1--;
822: }
823: m = k - m1;
824: m1 += 0x3FFF; // biased exponent of 2^(M1)
825: m += 0x3FFF; // biased exponent of 2^(M)
826:
827: goto expcont1;
828: }
829: }
830: } else { // |X| < 2^(-65)
831: float_rounding_mode = user_rnd_mode;
832: floatx80_rounding_precision = user_rnd_prec;
833:
834: a = floatx80_add(a, float32_to_floatx80(0x3F800000)); // 1 + X
835:
836: float_raise(float_flag_inexact);
837:
838: return a;
839: }
840: }
841:
842: /*----------------------------------------------------------------------------
843: | e to x minus 1
844: *----------------------------------------------------------------------------*/
845:
846: floatx80 floatx80_etoxm1(floatx80 a)
847: {
848: flag aSign;
849: int32 aExp;
850: bits64 aSig;
851:
852: int8 user_rnd_mode, user_rnd_prec;
853:
854: int32 compact, n, j, m, m1;
855: floatx80 fp0, fp1, fp2, fp3, l2, sc, onebysc;
856:
857: aSig = extractFloatx80Frac(a);
858: aExp = extractFloatx80Exp(a);
859: aSign = extractFloatx80Sign(a);
860:
861: if (aExp == 0x7FFF) {
862: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a);
863: if (aSign) return packFloatx80(aSign, one_exp, one_sig);
864: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low);
865: }
866:
867: if (aExp == 0 && aSig == 0) {
868: return packFloatx80(aSign, 0, 0);
869: }
870:
871: user_rnd_mode = float_rounding_mode;
872: user_rnd_prec = floatx80_rounding_precision;
873: float_rounding_mode = float_round_nearest_even;
874: floatx80_rounding_precision = 80;
875:
876: if (aExp >= 0x3FFD) { // |X| >= 1/4
877: compact = floatx80_make_compact(aExp, aSig);
878:
879: if (compact <= 0x4004C215) { // |X| <= 70 log2
880: fp0 = a;
881: fp1 = a;
882: fp0 = floatx80_mul(fp0, float32_to_floatx80(0x42B8AA3B)); // 64/log2 * X
883: n = floatx80_to_int32(fp0); // int(64/log2*X)
884: fp0 = int32_to_floatx80(n);
885:
886: j = n & 0x3F; // J = N mod 64
887: m = n / 64; // NOTE: this is really arithmetic right shift by 6
888: if (n < 0 && j) { // arithmetic right shift is division and round towards minus infinity
889: m--;
890: }
891: m1 = -m;
892: //m += 0x3FFF; // biased exponent of 2^(M)
893: //m1 += 0x3FFF; // biased exponent of -2^(-M)
894:
895: fp2 = fp0; // N
896: fp0 = floatx80_mul(fp0, float32_to_floatx80(0xBC317218)); // N * L1, L1 = lead(-log2/64)
897: l2 = packFloatx80(0, 0x3FDC, LIT64(0x82E308654361C4C6));
898: fp2 = floatx80_mul(fp2, l2); // N * L2, L1+L2 = -log2/64
899: fp0 = floatx80_add(fp0, fp1); // X + N*L1
900: fp0 = floatx80_add(fp0, fp2); // R
901:
902: fp1 = floatx80_mul(fp0, fp0); // S = R*R
903: fp2 = float32_to_floatx80(0x3950097B); // A6
904: fp2 = floatx80_mul(fp2, fp1); // fp2 is S*A6
905: fp3 = floatx80_mul(float32_to_floatx80(0x3AB60B6A), fp1); // fp3 is S*A5
906: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3F81111111174385))); // fp2 IS A4+S*A6
907: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0x3FA5555555554F5A))); // fp3 is A3+S*A5
908: fp2 = floatx80_mul(fp2, fp1); // fp2 IS S*(A4+S*A6)
909: fp3 = floatx80_mul(fp3, fp1); // fp3 IS S*(A3+S*A5)
910: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FC5555555555555))); // fp2 IS A2+S*(A4+S*A6)
911: fp3 = floatx80_add(fp3, float32_to_floatx80(0x3F000000)); // fp3 IS A1+S*(A3+S*A5)
912: fp2 = floatx80_mul(fp2, fp1); // fp2 IS S*(A2+S*(A4+S*A6))
913: fp1 = floatx80_mul(fp1, fp3); // fp1 IS S*(A1+S*(A3+S*A5))
914: fp2 = floatx80_mul(fp2, fp0); // fp2 IS R*S*(A2+S*(A4+S*A6))
915: fp0 = floatx80_add(fp0, fp1); // fp0 IS R+S*(A1+S*(A3+S*A5))
916: fp0 = floatx80_add(fp0, fp2); // fp0 IS EXP(R) - 1
917:
918: fp0 = floatx80_mul(fp0, exp_tbl[j]); // 2^(J/64)*(Exp(R)-1)
919:
920: if (m >= 64) {
921: fp1 = float32_to_floatx80(exp_tbl2[j]);
922: onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); // -2^(-M)
923: fp1 = floatx80_add(fp1, onebysc);
924: fp0 = floatx80_add(fp0, fp1);
925: fp0 = floatx80_add(fp0, exp_tbl[j]);
926: } else if (m < -3) {
927: fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j]));
928: fp0 = floatx80_add(fp0, exp_tbl[j]);
929: onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); // -2^(-M)
930: fp0 = floatx80_add(fp0, onebysc);
931: } else { // -3 <= m <= 63
932: fp1 = exp_tbl[j];
933: fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j]));
934: onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); // -2^(-M)
935: fp1 = floatx80_add(fp1, onebysc);
936: fp0 = floatx80_add(fp0, fp1);
937: }
938:
939: sc = packFloatx80(0, m + 0x3FFF, one_sig);
940:
941: float_rounding_mode = user_rnd_mode;
942: floatx80_rounding_precision = user_rnd_prec;
943:
944: a = floatx80_mul(fp0, sc);
945:
946: float_raise(float_flag_inexact);
947:
948: return a;
949: } else { // |X| > 70 log2
950: if (aSign) {
951: fp0 = float32_to_floatx80(0xBF800000); // -1
952:
953: float_rounding_mode = user_rnd_mode;
954: floatx80_rounding_precision = user_rnd_prec;
955:
956: a = floatx80_add(fp0, float32_to_floatx80(0x00800000)); // -1 + 2^(-126)
957:
958: float_raise(float_flag_inexact);
959:
960: return a;
961: } else {
962: float_rounding_mode = user_rnd_mode;
963: floatx80_rounding_precision = user_rnd_prec;
964:
965: return floatx80_etox(a);
966: }
967: }
968: } else { // |X| < 1/4
969: if (aExp >= 0x3FBE) {
970: fp0 = a;
971: fp0 = floatx80_mul(fp0, fp0); // S = X*X
972: fp1 = float32_to_floatx80(0x2F30CAA8); // B12
973: fp1 = floatx80_mul(fp1, fp0); // S * B12
974: fp2 = float32_to_floatx80(0x310F8290); // B11
975: fp1 = floatx80_add(fp1, float32_to_floatx80(0x32D73220)); // B10
976: fp2 = floatx80_mul(fp2, fp0);
977: fp1 = floatx80_mul(fp1, fp0);
978: fp2 = floatx80_add(fp2, float32_to_floatx80(0x3493F281)); // B9
979: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3EC71DE3A5774682))); // B8
980: fp2 = floatx80_mul(fp2, fp0);
981: fp1 = floatx80_mul(fp1, fp0);
982: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3EFA01A019D7CB68))); // B7
983: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3F2A01A01A019DF3))); // B6
984: fp2 = floatx80_mul(fp2, fp0);
985: fp1 = floatx80_mul(fp1, fp0);
986: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3F56C16C16C170E2))); // B5
987: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3F81111111111111))); // B4
988: fp2 = floatx80_mul(fp2, fp0);
989: fp1 = floatx80_mul(fp1, fp0);
990: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FA5555555555555))); // B3
991: fp3 = packFloatx80(0, 0x3FFC, LIT64(0xAAAAAAAAAAAAAAAB));
992: fp1 = floatx80_add(fp1, fp3); // B2
993: fp2 = floatx80_mul(fp2, fp0);
994: fp1 = floatx80_mul(fp1, fp0);
995:
996: fp2 = floatx80_mul(fp2, fp0);
997: fp1 = floatx80_mul(fp1, a);
998:
999: fp0 = floatx80_mul(fp0, float32_to_floatx80(0x3F000000)); // S*B1
1000: fp1 = floatx80_add(fp1, fp2); // Q
1001: fp0 = floatx80_add(fp0, fp1); // S*B1+Q
1002:
1003: float_rounding_mode = user_rnd_mode;
1004: floatx80_rounding_precision = user_rnd_prec;
1005:
1006: a = floatx80_add(fp0, a);
1007:
1008: float_raise(float_flag_inexact);
1009:
1010: return a;
1011: } else { // |X| < 2^(-65)
1012: sc = packFloatx80(1, 1, one_sig);
1013: fp0 = a;
1014:
1015: if (aExp < 0x0033) { // |X| < 2^(-16382)
1016: fp0 = floatx80_mul(fp0, float64_to_floatx80(LIT64(0x48B0000000000000)));
1017: fp0 = floatx80_add(fp0, sc);
1018:
1019: float_rounding_mode = user_rnd_mode;
1020: floatx80_rounding_precision = user_rnd_prec;
1021:
1022: a = floatx80_mul(fp0, float64_to_floatx80(LIT64(0x3730000000000000)));
1023: } else {
1024: float_rounding_mode = user_rnd_mode;
1025: floatx80_rounding_precision = user_rnd_prec;
1026:
1027: a = floatx80_add(fp0, sc);
1028: }
1029:
1030: float_raise(float_flag_inexact);
1031:
1032: return a;
1033: }
1034: }
1035: }
1036:
1037: /*----------------------------------------------------------------------------
1038: | Log base 10
1039: *----------------------------------------------------------------------------*/
1040:
1041: floatx80 floatx80_log10(floatx80 a)
1042: {
1043: flag aSign;
1044: int32 aExp;
1045: bits64 aSig;
1046:
1047: int8 user_rnd_mode, user_rnd_prec;
1048:
1049: floatx80 fp0, fp1;
1050:
1051: aSig = extractFloatx80Frac(a);
1052: aExp = extractFloatx80Exp(a);
1053: aSign = extractFloatx80Sign(a);
1054:
1055: if (aExp == 0x7FFF) {
1056: if ((bits64) (aSig<<1)) propagateFloatx80NaNOneArg(a);
1057: if (aSign == 0)
1058: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low);
1059: }
1060:
1061: if (aExp == 0 && aSig == 0) {
1062: float_raise(float_flag_divbyzero);
1063: return packFloatx80(1, 0x7FFF, floatx80_default_infinity_low);
1064: }
1065:
1066: if (aSign) {
1067: float_raise(float_flag_invalid);
1068: a.low = floatx80_default_nan_low;
1069: a.high = floatx80_default_nan_high;
1070: return a;
1071: }
1072:
1073: user_rnd_mode = float_rounding_mode;
1074: user_rnd_prec = floatx80_rounding_precision;
1075: float_rounding_mode = float_round_nearest_even;
1076: floatx80_rounding_precision = 80;
1077:
1078: fp0 = floatx80_logn(a);
1079: fp1 = packFloatx80(0, 0x3FFD, LIT64(0xDE5BD8A937287195)); // INV_L10
1080:
1081: float_rounding_mode = user_rnd_mode;
1082: floatx80_rounding_precision = user_rnd_prec;
1083:
1084: a = floatx80_mul(fp0, fp1); // LOGN(X)*INV_L10
1085:
1086: float_raise(float_flag_inexact);
1087:
1088: return a;
1089: }
1090:
1091: /*----------------------------------------------------------------------------
1092: | Log base 2
1093: *----------------------------------------------------------------------------*/
1094:
1095: floatx80 floatx80_log2(floatx80 a)
1096: {
1097: flag aSign;
1098: int32 aExp;
1099: bits64 aSig;
1100:
1101: int8 user_rnd_mode, user_rnd_prec;
1102:
1103: floatx80 fp0, fp1;
1104:
1105: aSig = extractFloatx80Frac(a);
1106: aExp = extractFloatx80Exp(a);
1107: aSign = extractFloatx80Sign(a);
1108:
1109: if (aExp == 0x7FFF) {
1110: if ((bits64) (aSig<<1)) propagateFloatx80NaNOneArg(a);
1111: if (aSign == 0)
1112: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low);
1113: }
1114:
1115: if (aExp == 0) {
1116: if (aSig == 0) {
1117: float_raise(float_flag_divbyzero);
1118: return packFloatx80(1, 0x7FFF, floatx80_default_infinity_low);
1119: }
1120: normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
1121: }
1122:
1123: if (aSign) {
1124: float_raise(float_flag_invalid);
1125: a.low = floatx80_default_nan_low;
1126: a.high = floatx80_default_nan_high;
1127: return a;
1128: }
1129:
1130: user_rnd_mode = float_rounding_mode;
1131: user_rnd_prec = floatx80_rounding_precision;
1132: float_rounding_mode = float_round_nearest_even;
1133: floatx80_rounding_precision = 80;
1134:
1135: if (aSig == one_sig) { // X is 2^k
1136: float_rounding_mode = user_rnd_mode;
1137: floatx80_rounding_precision = user_rnd_prec;
1138:
1139: a = int32_to_floatx80(aExp-0x3FFF);
1140: } else {
1141: fp0 = floatx80_logn(a);
1142: fp1 = packFloatx80(0, 0x3FFF, LIT64(0xB8AA3B295C17F0BC)); // INV_L2
1143:
1144: float_rounding_mode = user_rnd_mode;
1145: floatx80_rounding_precision = user_rnd_prec;
1146:
1147: a = floatx80_mul(fp0, fp1); // LOGN(X)*INV_L2
1148: }
1149:
1150: float_raise(float_flag_inexact);
1151:
1152: return a;
1153: }
1154:
1155: /*----------------------------------------------------------------------------
1156: | Log base e
1157: *----------------------------------------------------------------------------*/
1158:
1159: floatx80 floatx80_logn(floatx80 a)
1160: {
1161: flag aSign;
1162: int32 aExp;
1163: bits64 aSig, fSig;
1164:
1165: int8 user_rnd_mode, user_rnd_prec;
1166:
1167: int32 compact, j, k, adjk;
1168: floatx80 fp0, fp1, fp2, fp3, f, logof2, klog2, saveu;
1169:
1170: aSig = extractFloatx80Frac(a);
1171: aExp = extractFloatx80Exp(a);
1172: aSign = extractFloatx80Sign(a);
1173:
1174: if (aExp == 0x7FFF) {
1175: if ((bits64) (aSig<<1)) propagateFloatx80NaNOneArg(a);
1176: if (aSign == 0)
1177: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low);
1178: }
1179:
1180: adjk = 0;
1181:
1182: if (aExp == 0) {
1183: if (aSig == 0) { // zero
1184: float_raise(float_flag_divbyzero);
1185: return packFloatx80(1, 0x7FFF, floatx80_default_infinity_low);
1186: }
1187: #if 1
1188: if ((aSig & one_sig) == 0) { // denormal
1189: normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
1190: adjk = -100;
1191: aExp += 100;
1192: a = packFloatx80(aSign, aExp, aSig);
1193: }
1194: #else
1195: normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
1196: #endif
1197: }
1198:
1199: if (aSign) {
1200: float_raise(float_flag_invalid);
1201: a.low = floatx80_default_nan_low;
1202: a.high = floatx80_default_nan_high;
1203: return a;
1204: }
1205:
1206: user_rnd_mode = float_rounding_mode;
1207: user_rnd_prec = floatx80_rounding_precision;
1208: float_rounding_mode = float_round_nearest_even;
1209: floatx80_rounding_precision = 80;
1210:
1211: compact = floatx80_make_compact(aExp, aSig);
1212:
1213: if (compact < 0x3FFEF07D || compact > 0x3FFF8841) { // |X| < 15/16 or |X| > 17/16
1214: k = aExp - 0x3FFF;
1215: k += adjk;
1216: fp1 = int32_to_floatx80(k);
1217:
1218: fSig = (aSig & LIT64(0xFE00000000000000)) | LIT64(0x0100000000000000);
1219: j = (fSig >> 56) & 0x7E; // DISPLACEMENT FOR 1/F
1220:
1221: f = packFloatx80(0, 0x3FFF, fSig); // F
1222: fp0 = packFloatx80(0, 0x3FFF, aSig); // Y
1223:
1224: fp0 = floatx80_sub(fp0, f); // Y-F
1225:
1226: // LP1CONT1
1227: fp0 = floatx80_mul(fp0, log_tbl[j]); // FP0 IS U = (Y-F)/F
1228: logof2 = packFloatx80(0, 0x3FFE, LIT64(0xB17217F7D1CF79AC));
1229: klog2 = floatx80_mul(fp1, logof2); // FP1 IS K*LOG2
1230: fp2 = floatx80_mul(fp0, fp0); // FP2 IS V=U*U
1231:
1232: fp3 = fp2;
1233: fp1 = fp2;
1234:
1235: fp1 = floatx80_mul(fp1, float64_to_floatx80(LIT64(0x3FC2499AB5E4040B))); // V*A6
1236: fp2 = floatx80_mul(fp2, float64_to_floatx80(LIT64(0xBFC555B5848CB7DB))); // V*A5
1237: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3FC99999987D8730))); // A4+V*A6
1238: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0xBFCFFFFFFF6F7E97))); // A3+V*A5
1239: fp1 = floatx80_mul(fp1, fp3); // V*(A4+V*A6)
1240: fp2 = floatx80_mul(fp2, fp3); // V*(A3+V*A5)
1241: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3FD55555555555A4))); // A2+V*(A4+V*A6)
1242: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0xBFE0000000000008))); // A1+V*(A3+V*A5)
1243: fp1 = floatx80_mul(fp1, fp3); // V*(A2+V*(A4+V*A6))
1244: fp2 = floatx80_mul(fp2, fp3); // V*(A1+V*(A3+V*A5))
1245: fp1 = floatx80_mul(fp1, fp0); // U*V*(A2+V*(A4+V*A6))
1246: fp0 = floatx80_add(fp0, fp2); // U+V*(A1+V*(A3+V*A5))
1247:
1248: fp1 = floatx80_add(fp1, log_tbl[j+1]); // LOG(F)+U*V*(A2+V*(A4+V*A6))
1249: fp0 = floatx80_add(fp0, fp1); // FP0 IS LOG(F) + LOG(1+U)
1250:
1251: float_rounding_mode = user_rnd_mode;
1252: floatx80_rounding_precision = user_rnd_prec;
1253:
1254: a = floatx80_add(fp0, klog2);
1255:
1256: float_raise(float_flag_inexact);
1257:
1258: return a;
1259: } else { // |X-1| >= 1/16
1260: fp0 = a;
1261: fp1 = a;
1262: fp1 = floatx80_sub(fp1, float32_to_floatx80(0x3F800000)); // FP1 IS X-1
1263: fp0 = floatx80_add(fp0, float32_to_floatx80(0x3F800000)); // FP0 IS X+1
1264: fp1 = floatx80_add(fp1, fp1); // FP1 IS 2(X-1)
1265:
1266: // LP1CONT2
1267: fp1 = floatx80_div(fp1, fp0); // U
1268: saveu = fp1;
1269: fp0 = floatx80_mul(fp1, fp1); // FP0 IS V = U*U
1270: fp1 = floatx80_mul(fp0, fp0); // FP1 IS W = V*V
1271:
1272: fp3 = float64_to_floatx80(LIT64(0x3F175496ADD7DAD6)); // B5
1273: fp2 = float64_to_floatx80(LIT64(0x3F3C71C2FE80C7E0)); // B4
1274: fp3 = floatx80_mul(fp3, fp1); // W*B5
1275: fp2 = floatx80_mul(fp2, fp1); // W*B4
1276: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0x3F624924928BCCFF))); // B3+W*B5
1277: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3F899999999995EC))); // B2+W*B4
1278: fp1 = floatx80_mul(fp1, fp3); // W*(B3+W*B5)
1279: fp2 = floatx80_mul(fp2, fp0); // V*(B2+W*B4)
1280: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3FB5555555555555))); // B1+W*(B3+W*B5)
1281:
1282: fp0 = floatx80_mul(fp0, saveu); // FP0 IS U*V
1283: fp1 = floatx80_add(fp1, fp2); // B1+W*(B3+W*B5) + V*(B2+W*B4)
1284: fp0 = floatx80_mul(fp0, fp1); // U*V*( [B1+W*(B3+W*B5)] + [V*(B2+W*B4)] )
1285:
1286: float_rounding_mode = user_rnd_mode;
1287: floatx80_rounding_precision = user_rnd_prec;
1288:
1289: a = floatx80_add(fp0, saveu);
1290:
1291: //if (!floatx80_is_zero(a)) {
1292: float_raise(float_flag_inexact);
1293: //}
1294:
1295: return a;
1296: }
1297: }
1298:
1299: /*----------------------------------------------------------------------------
1300: | Log base e of x plus 1
1301: *----------------------------------------------------------------------------*/
1302:
1303: floatx80 floatx80_lognp1(floatx80 a)
1304: {
1305: flag aSign;
1306: int32 aExp;
1307: bits64 aSig, fSig;
1308:
1309: int8 user_rnd_mode, user_rnd_prec;
1310:
1311: int32 compact, j, k;
1312: floatx80 fp0, fp1, fp2, fp3, f, logof2, klog2, saveu;
1313:
1314: aSig = extractFloatx80Frac(a);
1315: aExp = extractFloatx80Exp(a);
1316: aSign = extractFloatx80Sign(a);
1317:
1318: if (aExp == 0x7FFF) {
1319: if ((bits64) (aSig<<1)) propagateFloatx80NaNOneArg(a);
1320: if (aSign) {
1321: float_raise(float_flag_invalid);
1322: a.low = floatx80_default_nan_low;
1323: a.high = floatx80_default_nan_high;
1324: return a;
1325: }
1326: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low);
1327: }
1328:
1329: if (aExp == 0 && aSig == 0) {
1330: return packFloatx80(aSign, 0, 0);
1331: }
1332:
1333: if (aSign && aExp >= one_exp) {
1334: if (aExp == one_exp && aSig == one_sig) {
1335: float_raise(float_flag_divbyzero);
1336: packFloatx80(aSign, 0x7FFF, floatx80_default_infinity_low);
1337: }
1338: float_raise(float_flag_invalid);
1339: a.low = floatx80_default_nan_low;
1340: a.high = floatx80_default_nan_high;
1341: return a;
1342: }
1343:
1344: if (aExp < 0x3f99 || (aExp == 0x3f99 && aSig == one_sig)) { // <= min threshold
1345: float_raise(float_flag_inexact);
1346: return floatx80_move(a);
1347: }
1348:
1349: user_rnd_mode = float_rounding_mode;
1350: user_rnd_prec = floatx80_rounding_precision;
1351: float_rounding_mode = float_round_nearest_even;
1352: floatx80_rounding_precision = 80;
1353:
1354: compact = floatx80_make_compact(aExp, aSig);
1355:
1356: fp0 = a; // Z
1357: fp1 = a;
1358:
1359: fp0 = floatx80_add(fp0, float32_to_floatx80(0x3F800000)); // X = (1+Z)
1360:
1361: aExp = extractFloatx80Exp(fp0);
1362: aSig = extractFloatx80Frac(fp0);
1363:
1364: compact = floatx80_make_compact(aExp, aSig);
1365:
1366: if (compact < 0x3FFE8000 || compact > 0x3FFFC000) { // |X| < 1/2 or |X| > 3/2
1367: k = aExp - 0x3FFF;
1368: fp1 = int32_to_floatx80(k);
1369:
1370: fSig = (aSig & LIT64(0xFE00000000000000)) | LIT64(0x0100000000000000);
1371: j = (fSig >> 56) & 0x7E; // DISPLACEMENT FOR 1/F
1372:
1373: f = packFloatx80(0, 0x3FFF, fSig); // F
1374: fp0 = packFloatx80(0, 0x3FFF, aSig); // Y
1375:
1376: fp0 = floatx80_sub(fp0, f); // Y-F
1377:
1378: lp1cont1:
1379: // LP1CONT1
1380: fp0 = floatx80_mul(fp0, log_tbl[j]); // FP0 IS U = (Y-F)/F
1381: logof2 = packFloatx80(0, 0x3FFE, LIT64(0xB17217F7D1CF79AC));
1382: klog2 = floatx80_mul(fp1, logof2); // FP1 IS K*LOG2
1383: fp2 = floatx80_mul(fp0, fp0); // FP2 IS V=U*U
1384:
1385: fp3 = fp2;
1386: fp1 = fp2;
1387:
1388: fp1 = floatx80_mul(fp1, float64_to_floatx80(LIT64(0x3FC2499AB5E4040B))); // V*A6
1389: fp2 = floatx80_mul(fp2, float64_to_floatx80(LIT64(0xBFC555B5848CB7DB))); // V*A5
1390: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3FC99999987D8730))); // A4+V*A6
1391: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0xBFCFFFFFFF6F7E97))); // A3+V*A5
1392: fp1 = floatx80_mul(fp1, fp3); // V*(A4+V*A6)
1393: fp2 = floatx80_mul(fp2, fp3); // V*(A3+V*A5)
1394: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3FD55555555555A4))); // A2+V*(A4+V*A6)
1395: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0xBFE0000000000008))); // A1+V*(A3+V*A5)
1396: fp1 = floatx80_mul(fp1, fp3); // V*(A2+V*(A4+V*A6))
1397: fp2 = floatx80_mul(fp2, fp3); // V*(A1+V*(A3+V*A5))
1398: fp1 = floatx80_mul(fp1, fp0); // U*V*(A2+V*(A4+V*A6))
1399: fp0 = floatx80_add(fp0, fp2); // U+V*(A1+V*(A3+V*A5))
1400:
1401: fp1 = floatx80_add(fp1, log_tbl[j+1]); // LOG(F)+U*V*(A2+V*(A4+V*A6))
1402: fp0 = floatx80_add(fp0, fp1); // FP0 IS LOG(F) + LOG(1+U)
1403:
1404: float_rounding_mode = user_rnd_mode;
1405: floatx80_rounding_precision = user_rnd_prec;
1406:
1407: a = floatx80_add(fp0, klog2);
1408:
1409: float_raise(float_flag_inexact);
1410:
1411: return a;
1412: } else if (compact < 0x3FFEF07D || compact > 0x3FFF8841) { // |X| < 1/16 or |X| > -1/16
1413: // LP1CARE
1414: fSig = (aSig & LIT64(0xFE00000000000000)) | LIT64(0x0100000000000000);
1415: f = packFloatx80(0, 0x3FFF, fSig); // F
1416: j = (fSig >> 56) & 0x7E; // DISPLACEMENT FOR 1/F
1417:
1418: if (compact >= 0x3FFF8000) { // 1+Z >= 1
1419: // KISZERO
1420: fp0 = floatx80_sub(float32_to_floatx80(0x3F800000), f); // 1-F
1421: fp0 = floatx80_add(fp0, fp1); // FP0 IS Y-F = (1-F)+Z
1422: fp1 = packFloatx80(0, 0, 0); // K = 0
1423: } else {
1424: // KISNEG
1425: fp0 = floatx80_sub(float32_to_floatx80(0x40000000), f); // 2-F
1426: fp1 = floatx80_add(fp1, fp1); // 2Z
1427: fp0 = floatx80_add(fp0, fp1); // FP0 IS Y-F = (2-F)+2Z
1428: fp1 = packFloatx80(1, one_exp, one_sig); // K = -1
1429: }
1430: goto lp1cont1;
1431: } else {
1432: // LP1ONE16
1433: fp1 = floatx80_add(fp1, fp1); // FP1 IS 2Z
1434: fp0 = floatx80_add(fp0, float32_to_floatx80(0x3F800000)); // FP0 IS 1+X
1435:
1436: // LP1CONT2
1437: fp1 = floatx80_div(fp1, fp0); // U
1438: saveu = fp1;
1439: fp0 = floatx80_mul(fp1, fp1); // FP0 IS V = U*U
1440: fp1 = floatx80_mul(fp0, fp0); // FP1 IS W = V*V
1441:
1442: fp3 = float64_to_floatx80(LIT64(0x3F175496ADD7DAD6)); // B5
1443: fp2 = float64_to_floatx80(LIT64(0x3F3C71C2FE80C7E0)); // B4
1444: fp3 = floatx80_mul(fp3, fp1); // W*B5
1445: fp2 = floatx80_mul(fp2, fp1); // W*B4
1446: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0x3F624924928BCCFF))); // B3+W*B5
1447: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3F899999999995EC))); // B2+W*B4
1448: fp1 = floatx80_mul(fp1, fp3); // W*(B3+W*B5)
1449: fp2 = floatx80_mul(fp2, fp0); // V*(B2+W*B4)
1450: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3FB5555555555555))); // B1+W*(B3+W*B5)
1451:
1452: fp0 = floatx80_mul(fp0, saveu); // FP0 IS U*V
1453: fp1 = floatx80_add(fp1, fp2); // B1+W*(B3+W*B5) + V*(B2+W*B4)
1454: fp0 = floatx80_mul(fp0, fp1); // U*V*( [B1+W*(B3+W*B5)] + [V*(B2+W*B4)] )
1455:
1456: float_rounding_mode = user_rnd_mode;
1457: floatx80_rounding_precision = user_rnd_prec;
1458:
1459: a = floatx80_add(fp0, saveu);
1460:
1461: //if (!floatx80_is_zero(a)) {
1462: float_raise(float_flag_inexact);
1463: //}
1464:
1465: return a;
1466: }
1467: }
1468:
1469: /*----------------------------------------------------------------------------
1470: | Sine
1471: *----------------------------------------------------------------------------*/
1472:
1473: floatx80 floatx80_sin(floatx80 a)
1474: {
1475: flag aSign, xSign;
1476: int32 aExp, xExp;
1477: bits64 aSig, xSig;
1478:
1479: int8 user_rnd_mode, user_rnd_prec;
1480:
1481: int32 compact, l, n, j;
1482: floatx80 fp0, fp1, fp2, fp3, fp4, fp5, x, invtwopi, twopi1, twopi2;
1483: float32 posneg1, twoto63;
1484: flag adjn, endflag;
1485:
1486: aSig = extractFloatx80Frac(a);
1487: aExp = extractFloatx80Exp(a);
1488: aSign = extractFloatx80Sign(a);
1489:
1490: if (aExp == 0x7FFF) {
1491: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a);
1492: float_raise(float_flag_invalid);
1493: a.low = floatx80_default_nan_low;
1494: a.high = floatx80_default_nan_high;
1495: return a;
1496: }
1497:
1498: if (aExp == 0 && aSig == 0) {
1499: return packFloatx80(aSign, 0, 0);
1500: }
1501:
1502: adjn = 0;
1503:
1504: user_rnd_mode = float_rounding_mode;
1505: user_rnd_prec = floatx80_rounding_precision;
1506: float_rounding_mode = float_round_nearest_even;
1507: floatx80_rounding_precision = 80;
1508:
1509: compact = floatx80_make_compact(aExp, aSig);
1510:
1511: fp0 = a;
1512:
1513: if (compact < 0x3FD78000 || compact > 0x4004BC7E) { // 2^(-40) > |X| > 15 PI
1514: if (compact > 0x3FFF8000) { // |X| >= 15 PI
1515: // REDUCEX
1516: fp1 = packFloatx80(0, 0, 0);
1517: if (compact == 0x7FFEFFFF) {
1518: twopi1 = packFloatx80(aSign ^ 1, 0x7FFE, LIT64(0xC90FDAA200000000));
1519: twopi2 = packFloatx80(aSign ^ 1, 0x7FDC, LIT64(0x85A308D300000000));
1520: fp0 = floatx80_add(fp0, twopi1);
1521: fp1 = fp0;
1522: fp0 = floatx80_add(fp0, twopi2);
1523: fp1 = floatx80_sub(fp1, fp0);
1524: fp1 = floatx80_add(fp1, twopi2);
1525: }
1526: loop:
1527: xSign = extractFloatx80Sign(fp0);
1528: xExp = extractFloatx80Exp(fp0);
1529: xExp -= 0x3FFF;
1530: if (xExp <= 28) {
1531: l = 0;
1532: endflag = 1;
1533: } else {
1534: l = xExp - 27;
1535: endflag = 0;
1536: }
1537: invtwopi = packFloatx80(0, 0x3FFE - l, LIT64(0xA2F9836E4E44152A)); // INVTWOPI
1538: twopi1 = packFloatx80(0, 0x3FFF + l, LIT64(0xC90FDAA200000000));
1539: twopi2 = packFloatx80(0, 0x3FDD + l, LIT64(0x85A308D300000000));
1540:
1541: twoto63 = 0x5F000000;
1542: twoto63 |= xSign ? 0x80000000 : 0x00000000; // SIGN(INARG)*2^63 IN SGL
1543:
1544: fp2 = floatx80_mul(fp0, invtwopi);
1545: fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63)); // THE FRACTIONAL PART OF FP2 IS ROUNDED
1546: fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63)); // FP2 is N
1547: fp4 = floatx80_mul(twopi1, fp2); // W = N*P1
1548: fp5 = floatx80_mul(twopi2, fp2); // w = N*P2
1549: fp3 = floatx80_add(fp4, fp5); // FP3 is P
1550: fp4 = floatx80_sub(fp4, fp3); // W-P
1551: fp0 = floatx80_sub(fp0, fp3); // FP0 is A := R - P
1552: fp4 = floatx80_add(fp4, fp5); // FP4 is p = (W-P)+w
1553: fp3 = fp0; // FP3 is A
1554: fp1 = floatx80_sub(fp1, fp4); // FP1 is a := r - p
1555: fp0 = floatx80_add(fp0, fp1); // FP0 is R := A+a
1556:
1557: if (endflag > 0) {
1558: n = floatx80_to_int32(fp2);
1559: goto sincont;
1560: }
1561: fp3 = floatx80_sub(fp3, fp0); // A-R
1562: fp1 = floatx80_add(fp1, fp3); // FP1 is r := (A-R)+a
1563: goto loop;
1564: } else {
1565: // SINSM
1566: fp0 = float32_to_floatx80(0x3F800000); // 1
1567:
1568: float_rounding_mode = user_rnd_mode;
1569: floatx80_rounding_precision = user_rnd_prec;
1570:
1571: if (adjn) {
1572: // COSTINY
1573: a = floatx80_sub(fp0, float32_to_floatx80(0x00800000));
1574: } else {
1575: // SINTINY
1576: a = floatx80_move(a);
1577: }
1578: float_raise(float_flag_inexact);
1579:
1580: return a;
1581: }
1582: } else {
1583: fp1 = floatx80_mul(fp0, float64_to_floatx80(LIT64(0x3FE45F306DC9C883))); // X*2/PI
1584:
1585: n = floatx80_to_int32(fp1);
1586: j = 32 + n;
1587:
1588: fp0 = floatx80_sub(fp0, pi_tbl[j]); // X-Y1
1589: fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j])); // FP0 IS R = (X-Y1)-Y2
1590:
1591: sincont:
1592: if ((n + adjn) & 1) {
1593: // COSPOLY
1594: fp0 = floatx80_mul(fp0, fp0); // FP0 IS S
1595: fp1 = floatx80_mul(fp0, fp0); // FP1 IS T
1596: fp2 = float64_to_floatx80(LIT64(0x3D2AC4D0D6011EE3)); // B8
1597: fp3 = float64_to_floatx80(LIT64(0xBDA9396F9F45AC19)); // B7
1598:
1599: xSign = extractFloatx80Sign(fp0); // X IS S
1600: xExp = extractFloatx80Exp(fp0);
1601: xSig = extractFloatx80Frac(fp0);
1602:
1603: if (((n + adjn) >> 1) & 1) {
1604: xSign ^= 1;
1605: posneg1 = 0xBF800000; // -1
1606: } else {
1607: xSign ^= 0;
1608: posneg1 = 0x3F800000; // 1
1609: } // X IS NOW R'= SGN*R
1610:
1611: fp2 = floatx80_mul(fp2, fp1); // TB8
1612: fp3 = floatx80_mul(fp3, fp1); // TB7
1613: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3E21EED90612C972))); // B6+TB8
1614: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBE927E4FB79D9FCF))); // B5+TB7
1615: fp2 = floatx80_mul(fp2, fp1); // T(B6+TB8)
1616: fp3 = floatx80_mul(fp3, fp1); // T(B5+TB7)
1617: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3EFA01A01A01D423))); // B4+T(B6+TB8)
1618: fp4 = packFloatx80(1, 0x3FF5, LIT64(0xB60B60B60B61D438));
1619: fp3 = floatx80_add(fp3, fp4); // B3+T(B5+TB7)
1620: fp2 = floatx80_mul(fp2, fp1); // T(B4+T(B6+TB8))
1621: fp1 = floatx80_mul(fp1, fp3); // T(B3+T(B5+TB7))
1622: fp4 = packFloatx80(0, 0x3FFA, LIT64(0xAAAAAAAAAAAAAB5E));
1623: fp2 = floatx80_add(fp2, fp4); // B2+T(B4+T(B6+TB8))
1624: fp1 = floatx80_add(fp1, float32_to_floatx80(0xBF000000)); // B1+T(B3+T(B5+TB7))
1625: fp0 = floatx80_mul(fp0, fp2); // S(B2+T(B4+T(B6+TB8)))
1626: fp0 = floatx80_add(fp0, fp1); // [B1+T(B3+T(B5+TB7))]+[S(B2+T(B4+T(B6+TB8)))]
1627:
1628: x = packFloatx80(xSign, xExp, xSig);
1629: fp0 = floatx80_mul(fp0, x);
1630:
1631: float_rounding_mode = user_rnd_mode;
1632: floatx80_rounding_precision = user_rnd_prec;
1633:
1634: a = floatx80_add(fp0, float32_to_floatx80(posneg1));
1635:
1636: float_raise(float_flag_inexact);
1637:
1638: return a;
1639: } else {
1640: // SINPOLY
1641: xSign = extractFloatx80Sign(fp0); // X IS R
1642: xExp = extractFloatx80Exp(fp0);
1643: xSig = extractFloatx80Frac(fp0);
1644:
1645: xSign ^= ((n + adjn) >> 1) & 1; // X IS NOW R'= SGN*R
1646:
1647: fp0 = floatx80_mul(fp0, fp0); // FP0 IS S
1648: fp1 = floatx80_mul(fp0, fp0); // FP1 IS T
1649: fp3 = float64_to_floatx80(LIT64(0xBD6AAA77CCC994F5)); // A7
1650: fp2 = float64_to_floatx80(LIT64(0x3DE612097AAE8DA1)); // A6
1651: fp3 = floatx80_mul(fp3, fp1); // T*A7
1652: fp2 = floatx80_mul(fp2, fp1); // T*A6
1653: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBE5AE6452A118AE4))); // A5+T*A7
1654: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3EC71DE3A5341531))); // A4+T*A6
1655: fp3 = floatx80_mul(fp3, fp1); // T(A5+TA7)
1656: fp2 = floatx80_mul(fp2, fp1); // T(A4+TA6)
1657: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBF2A01A01A018B59))); // A3+T(A5+TA7)
1658: fp4 = packFloatx80(0, 0x3FF8, LIT64(0x88888888888859AF));
1659: fp2 = floatx80_add(fp2, fp4); // A2+T(A4+TA6)
1660: fp1 = floatx80_mul(fp1, fp3); // T(A3+T(A5+TA7))
1661: fp2 = floatx80_mul(fp2, fp0); // S(A2+T(A4+TA6))
1662: fp4 = packFloatx80(1, 0x3FFC, LIT64(0xAAAAAAAAAAAAAA99));
1663: fp1 = floatx80_add(fp1, fp4); // A1+T(A3+T(A5+TA7))
1664: fp1 = floatx80_add(fp1, fp2); // [A1+T(A3+T(A5+TA7))]+[S(A2+T(A4+TA6))]
1665:
1666: x = packFloatx80(xSign, xExp, xSig);
1667: fp0 = floatx80_mul(fp0, x); // R'*S
1668: fp0 = floatx80_mul(fp0, fp1); // SIN(R')-R'
1669:
1670: float_rounding_mode = user_rnd_mode;
1671: floatx80_rounding_precision = user_rnd_prec;
1672:
1673: a = floatx80_add(fp0, x);
1674:
1675: float_raise(float_flag_inexact);
1676:
1677: return a;
1678: }
1679: }
1680: }
1681:
1682: /*----------------------------------------------------------------------------
1683: | Hyperbolic sine
1684: *----------------------------------------------------------------------------*/
1685:
1686: floatx80 floatx80_sinh(floatx80 a)
1687: {
1688: flag aSign;
1689: int32 aExp;
1690: bits64 aSig;
1691:
1692: int8 user_rnd_mode, user_rnd_prec;
1693:
1694: int32 compact;
1695: floatx80 fp0, fp1, fp2;
1696: float32 fact;
1697:
1698: aSig = extractFloatx80Frac(a);
1699: aExp = extractFloatx80Exp(a);
1700: aSign = extractFloatx80Sign(a);
1701:
1702: if (aExp == 0x7FFF) {
1703: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a);
1704: return packFloatx80(aSign, 0x7FFF, floatx80_default_infinity_low);
1705: }
1706:
1707: if (aExp == 0 && aSig == 0) {
1708: return packFloatx80(aSign, 0, 0);
1709: }
1710:
1711: user_rnd_mode = float_rounding_mode;
1712: user_rnd_prec = floatx80_rounding_precision;
1713: float_rounding_mode = float_round_nearest_even;
1714: floatx80_rounding_precision = 80;
1715:
1716: compact = floatx80_make_compact(aExp, aSig);
1717:
1718: if (compact > 0x400CB167) {
1719: // SINHBIG
1720: if (compact > 0x400CB2B3) {
1721: float_rounding_mode = user_rnd_mode;
1722: floatx80_rounding_precision = user_rnd_prec;
1723:
1724: return roundAndPackFloatx80(floatx80_rounding_precision, aSign, 0x8000, aSig, 0);
1725: } else {
1726: fp0 = floatx80_abs(a); // Y = |X|
1727: fp0 = floatx80_sub(fp0, float64_to_floatx80(LIT64(0x40C62D38D3D64634))); // (|X|-16381LOG2_LEAD)
1728: fp0 = floatx80_sub(fp0, float64_to_floatx80(LIT64(0x3D6F90AEB1E75CC7))); // |X| - 16381 LOG2, ACCURATE
1729: fp0 = floatx80_etox(fp0);
1730: fp2 = packFloatx80(aSign, 0x7FFB, one_sig);
1731:
1732: float_rounding_mode = user_rnd_mode;
1733: floatx80_rounding_precision = user_rnd_prec;
1734:
1735: a = floatx80_mul(fp0, fp2);
1736:
1737: float_raise(float_flag_inexact);
1738:
1739: return a;
1740: }
1741: } else { // |X| < 16380 LOG2
1742: fp0 = floatx80_abs(a); // Y = |X|
1743: fp0 = floatx80_etoxm1(fp0); // FP0 IS Z = EXPM1(Y)
1744: fp1 = floatx80_add(fp0, float32_to_floatx80(0x3F800000)); // 1+Z
1745: fp2 = fp0;
1746: fp0 = floatx80_div(fp0, fp1); // Z/(1+Z)
1747: fp0 = floatx80_add(fp0, fp2);
1748:
1749: fact = 0x3F000000;
1750: fact |= aSign ? 0x80000000 : 0x00000000;
1751:
1752: float_rounding_mode = user_rnd_mode;
1753: floatx80_rounding_precision = user_rnd_prec;
1754:
1755: a = floatx80_mul(fp0, float32_to_floatx80(fact));
1756:
1757: float_raise(float_flag_inexact);
1758:
1759: return a;
1760: }
1761: }
1762:
1763: /*----------------------------------------------------------------------------
1764: | Tangent
1765: *----------------------------------------------------------------------------*/
1766:
1767: floatx80 floatx80_tan(floatx80 a)
1768: {
1769: flag aSign, xSign;
1770: int32 aExp, xExp;
1771: bits64 aSig, xSig;
1772:
1773: int8 user_rnd_mode, user_rnd_prec;
1774:
1775: int32 compact, l, n, j;
1776: floatx80 fp0, fp1, fp2, fp3, fp4, fp5, invtwopi, twopi1, twopi2;
1777: float32 twoto63;
1778: flag endflag;
1779:
1780: aSig = extractFloatx80Frac(a);
1781: aExp = extractFloatx80Exp(a);
1782: aSign = extractFloatx80Sign(a);
1783:
1784: if (aExp == 0x7FFF) {
1785: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a);
1786: float_raise(float_flag_invalid);
1787: a.low = floatx80_default_nan_low;
1788: a.high = floatx80_default_nan_high;
1789: return a;
1790: }
1791:
1792: if (aExp == 0 && aSig == 0) {
1793: return packFloatx80(aSign, 0, 0);
1794: }
1795:
1796: user_rnd_mode = float_rounding_mode;
1797: user_rnd_prec = floatx80_rounding_precision;
1798: float_rounding_mode = float_round_nearest_even;
1799: floatx80_rounding_precision = 80;
1800:
1801: compact = floatx80_make_compact(aExp, aSig);
1802:
1803: fp0 = a;
1804:
1805: if (compact < 0x3FD78000 || compact > 0x4004BC7E) { // 2^(-40) > |X| > 15 PI
1806: if (compact > 0x3FFF8000) { // |X| >= 15 PI
1807: // REDUCEX
1808: fp1 = packFloatx80(0, 0, 0);
1809: if (compact == 0x7FFEFFFF) {
1810: twopi1 = packFloatx80(aSign ^ 1, 0x7FFE, LIT64(0xC90FDAA200000000));
1811: twopi2 = packFloatx80(aSign ^ 1, 0x7FDC, LIT64(0x85A308D300000000));
1812: fp0 = floatx80_add(fp0, twopi1);
1813: fp1 = fp0;
1814: fp0 = floatx80_add(fp0, twopi2);
1815: fp1 = floatx80_sub(fp1, fp0);
1816: fp1 = floatx80_add(fp1, twopi2);
1817: }
1818: loop:
1819: xSign = extractFloatx80Sign(fp0);
1820: xExp = extractFloatx80Exp(fp0);
1821: xExp -= 0x3FFF;
1822: if (xExp <= 28) {
1823: l = 0;
1824: endflag = 1;
1825: } else {
1826: l = xExp - 27;
1827: endflag = 0;
1828: }
1829: invtwopi = packFloatx80(0, 0x3FFE - l, LIT64(0xA2F9836E4E44152A)); // INVTWOPI
1830: twopi1 = packFloatx80(0, 0x3FFF + l, LIT64(0xC90FDAA200000000));
1831: twopi2 = packFloatx80(0, 0x3FDD + l, LIT64(0x85A308D300000000));
1832:
1833: twoto63 = 0x5F000000;
1834: twoto63 |= xSign ? 0x80000000 : 0x00000000; // SIGN(INARG)*2^63 IN SGL
1835:
1836: fp2 = floatx80_mul(fp0, invtwopi);
1837: fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63)); // THE FRACTIONAL PART OF FP2 IS ROUNDED
1838: fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63)); // FP2 is N
1839: fp4 = floatx80_mul(twopi1, fp2); // W = N*P1
1840: fp5 = floatx80_mul(twopi2, fp2); // w = N*P2
1841: fp3 = floatx80_add(fp4, fp5); // FP3 is P
1842: fp4 = floatx80_sub(fp4, fp3); // W-P
1843: fp0 = floatx80_sub(fp0, fp3); // FP0 is A := R - P
1844: fp4 = floatx80_add(fp4, fp5); // FP4 is p = (W-P)+w
1845: fp3 = fp0; // FP3 is A
1846: fp1 = floatx80_sub(fp1, fp4); // FP1 is a := r - p
1847: fp0 = floatx80_add(fp0, fp1); // FP0 is R := A+a
1848:
1849: if (endflag > 0) {
1850: n = floatx80_to_int32(fp2);
1851: goto tancont;
1852: }
1853: fp3 = floatx80_sub(fp3, fp0); // A-R
1854: fp1 = floatx80_add(fp1, fp3); // FP1 is r := (A-R)+a
1855: goto loop;
1856: } else {
1857: float_rounding_mode = user_rnd_mode;
1858: floatx80_rounding_precision = user_rnd_prec;
1859:
1860: a = floatx80_move(a);
1861:
1862: float_raise(float_flag_inexact);
1863:
1864: return a;
1865: }
1866: } else {
1867: fp1 = floatx80_mul(fp0, float64_to_floatx80(LIT64(0x3FE45F306DC9C883))); // X*2/PI
1868:
1869: n = floatx80_to_int32(fp1);
1870: j = 32 + n;
1871:
1872: fp0 = floatx80_sub(fp0, pi_tbl[j]); // X-Y1
1873: fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j])); // FP0 IS R = (X-Y1)-Y2
1874:
1875: tancont:
1876: if (n & 1) {
1877: // NODD
1878: fp1 = fp0; // R
1879: fp0 = floatx80_mul(fp0, fp0); // S = R*R
1880: fp3 = float64_to_floatx80(LIT64(0x3EA0B759F50F8688)); // Q4
1881: fp2 = float64_to_floatx80(LIT64(0xBEF2BAA5A8924F04)); // P3
1882: fp3 = floatx80_mul(fp3, fp0); // SQ4
1883: fp2 = floatx80_mul(fp2, fp0); // SP3
1884: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBF346F59B39BA65F))); // Q3+SQ4
1885: fp4 = packFloatx80(0, 0x3FF6, LIT64(0xE073D3FC199C4A00));
1886: fp2 = floatx80_add(fp2, fp4); // P2+SP3
1887: fp3 = floatx80_mul(fp3, fp0); // S(Q3+SQ4)
1888: fp2 = floatx80_mul(fp2, fp0); // S(P2+SP3)
1889: fp4 = packFloatx80(0, 0x3FF9, LIT64(0xD23CD68415D95FA1));
1890: fp3 = floatx80_add(fp3, fp4); // Q2+S(Q3+SQ4)
1891: fp4 = packFloatx80(1, 0x3FFC, LIT64(0x8895A6C5FB423BCA));
1892: fp2 = floatx80_add(fp2, fp4); // P1+S(P2+SP3)
1893: fp3 = floatx80_mul(fp3, fp0); // S(Q2+S(Q3+SQ4))
1894: fp2 = floatx80_mul(fp2, fp0); // S(P1+S(P2+SP3))
1895: fp4 = packFloatx80(1, 0x3FFD, LIT64(0xEEF57E0DA84BC8CE));
1896: fp3 = floatx80_add(fp3, fp4); // Q1+S(Q2+S(Q3+SQ4))
1897: fp2 = floatx80_mul(fp2, fp1); // RS(P1+S(P2+SP3))
1898: fp0 = floatx80_mul(fp0, fp3); // S(Q1+S(Q2+S(Q3+SQ4)))
1899: fp1 = floatx80_add(fp1, fp2); // R+RS(P1+S(P2+SP3))
1900: fp0 = floatx80_add(fp0, float32_to_floatx80(0x3F800000)); // 1+S(Q1+S(Q2+S(Q3+SQ4)))
1901:
1902: xSign = extractFloatx80Sign(fp1);
1903: xExp = extractFloatx80Exp(fp1);
1904: xSig = extractFloatx80Frac(fp1);
1905: xSign ^= 1;
1906: fp1 = packFloatx80(xSign, xExp, xSig);
1907:
1908: float_rounding_mode = user_rnd_mode;
1909: floatx80_rounding_precision = user_rnd_prec;
1910:
1911: a = floatx80_div(fp0, fp1);
1912:
1913: float_raise(float_flag_inexact);
1914:
1915: return a;
1916: } else {
1917: fp1 = floatx80_mul(fp0, fp0); // S = R*R
1918: fp3 = float64_to_floatx80(LIT64(0x3EA0B759F50F8688)); // Q4
1919: fp2 = float64_to_floatx80(LIT64(0xBEF2BAA5A8924F04)); // P3
1920: fp3 = floatx80_mul(fp3, fp1); // SQ4
1921: fp2 = floatx80_mul(fp2, fp1); // SP3
1922: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBF346F59B39BA65F))); // Q3+SQ4
1923: fp4 = packFloatx80(0, 0x3FF6, LIT64(0xE073D3FC199C4A00));
1924: fp2 = floatx80_add(fp2, fp4); // P2+SP3
1925: fp3 = floatx80_mul(fp3, fp1); // S(Q3+SQ4)
1926: fp2 = floatx80_mul(fp2, fp1); // S(P2+SP3)
1927: fp4 = packFloatx80(0, 0x3FF9, LIT64(0xD23CD68415D95FA1));
1928: fp3 = floatx80_add(fp3, fp4); // Q2+S(Q3+SQ4)
1929: fp4 = packFloatx80(1, 0x3FFC, LIT64(0x8895A6C5FB423BCA));
1930: fp2 = floatx80_add(fp2, fp4); // P1+S(P2+SP3)
1931: fp3 = floatx80_mul(fp3, fp1); // S(Q2+S(Q3+SQ4))
1932: fp2 = floatx80_mul(fp2, fp1); // S(P1+S(P2+SP3))
1933: fp4 = packFloatx80(1, 0x3FFD, LIT64(0xEEF57E0DA84BC8CE));
1934: fp3 = floatx80_add(fp3, fp4); // Q1+S(Q2+S(Q3+SQ4))
1935: fp2 = floatx80_mul(fp2, fp0); // RS(P1+S(P2+SP3))
1936: fp1 = floatx80_mul(fp1, fp3); // S(Q1+S(Q2+S(Q3+SQ4)))
1937: fp0 = floatx80_add(fp0, fp2); // R+RS(P1+S(P2+SP3))
1938: fp1 = floatx80_add(fp1, float32_to_floatx80(0x3F800000)); // 1+S(Q1+S(Q2+S(Q3+SQ4)))
1939:
1940: float_rounding_mode = user_rnd_mode;
1941: floatx80_rounding_precision = user_rnd_prec;
1942:
1943: a = floatx80_div(fp0, fp1);
1944:
1945: float_raise(float_flag_inexact);
1946:
1947: return a;
1948: }
1949: }
1950: }
1951:
1952: /*----------------------------------------------------------------------------
1953: | Hyperbolic tangent
1954: *----------------------------------------------------------------------------*/
1955:
1956: floatx80 floatx80_tanh(floatx80 a)
1957: {
1958: flag aSign, vSign;
1959: int32 aExp, vExp;
1960: bits64 aSig, vSig;
1961:
1962: int8 user_rnd_mode, user_rnd_prec;
1963:
1964: int32 compact;
1965: floatx80 fp0, fp1;
1966: float32 sign;
1967:
1968: aSig = extractFloatx80Frac(a);
1969: aExp = extractFloatx80Exp(a);
1970: aSign = extractFloatx80Sign(a);
1971:
1972: if (aExp == 0x7FFF) {
1973: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a);
1974: return packFloatx80(aSign, one_exp, one_sig);
1975: }
1976:
1977: if (aExp == 0 && aSig == 0) {
1978: return packFloatx80(aSign, 0, 0);
1979: }
1980:
1981: user_rnd_mode = float_rounding_mode;
1982: user_rnd_prec = floatx80_rounding_precision;
1983: float_rounding_mode = float_round_nearest_even;
1984: floatx80_rounding_precision = 80;
1985:
1986: compact = floatx80_make_compact(aExp, aSig);
1987:
1988: if (compact < 0x3FD78000 || compact > 0x3FFFDDCE) {
1989: // TANHBORS
1990: if (compact < 0x3FFF8000) {
1991: // TANHSM
1992: float_rounding_mode = user_rnd_mode;
1993: floatx80_rounding_precision = user_rnd_prec;
1994:
1995: a = floatx80_move(a);
1996:
1997: float_raise(float_flag_inexact);
1998:
1999: return a;
2000: } else {
2001: if (compact > 0x40048AA1) {
2002: // TANHHUGE
2003: sign = 0x3F800000;
2004: sign |= aSign ? 0x80000000 : 0x00000000;
2005: fp0 = float32_to_floatx80(sign);
2006: sign &= 0x80000000;
2007: sign ^= 0x80800000; // -SIGN(X)*EPS
2008:
2009: float_rounding_mode = user_rnd_mode;
2010: floatx80_rounding_precision = user_rnd_prec;
2011:
2012: a = floatx80_add(fp0, float32_to_floatx80(sign));
2013:
2014: float_raise(float_flag_inexact);
2015:
2016: return a;
2017: } else {
2018: fp0 = packFloatx80(0, aExp+1, aSig); // Y = 2|X|
2019: fp0 = floatx80_etox(fp0); // FP0 IS EXP(Y)
2020: fp0 = floatx80_add(fp0, float32_to_floatx80(0x3F800000)); // EXP(Y)+1
2021: sign = aSign ? 0x80000000 : 0x00000000;
2022: fp1 = floatx80_div(float32_to_floatx80(sign^0xC0000000), fp0); // -SIGN(X)*2 / [EXP(Y)+1]
2023: fp0 = float32_to_floatx80(sign | 0x3F800000); // SIGN
2024:
2025: float_rounding_mode = user_rnd_mode;
2026: floatx80_rounding_precision = user_rnd_prec;
2027:
2028: a = floatx80_add(fp1, fp0);
2029:
2030: float_raise(float_flag_inexact);
2031:
2032: return a;
2033: }
2034: }
2035: } else { // 2**(-40) < |X| < (5/2)LOG2
2036: fp0 = packFloatx80(0, aExp+1, aSig); // Y = 2|X|
2037: fp0 = floatx80_etoxm1(fp0); // FP0 IS Z = EXPM1(Y)
2038: fp1 = floatx80_add(fp0, float32_to_floatx80(0x40000000)); // Z+2
2039:
2040: vSign = extractFloatx80Sign(fp1);
2041: vExp = extractFloatx80Exp(fp1);
2042: vSig = extractFloatx80Frac(fp1);
2043:
2044: fp1 = packFloatx80(vSign ^ aSign, vExp, vSig);
2045:
2046: float_rounding_mode = user_rnd_mode;
2047: floatx80_rounding_precision = user_rnd_prec;
2048:
2049: a = floatx80_div(fp0, fp1);
2050:
2051: float_raise(float_flag_inexact);
2052:
2053: return a;
2054: }
2055: }
2056:
2057: /*----------------------------------------------------------------------------
2058: | 10 to x
2059: *----------------------------------------------------------------------------*/
2060:
2061: floatx80 floatx80_tentox(floatx80 a)
2062: {
2063: flag aSign;
2064: int32 aExp;
2065: bits64 aSig;
2066:
2067: int8 user_rnd_mode, user_rnd_prec;
2068:
2069: int32 compact, n, j, l, m, m1;
2070: floatx80 fp0, fp1, fp2, fp3, adjfact, fact1, fact2;
2071:
2072: aSig = extractFloatx80Frac(a);
2073: aExp = extractFloatx80Exp(a);
2074: aSign = extractFloatx80Sign(a);
2075:
2076: if (aExp == 0x7FFF) {
2077: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a);
2078: if (aSign) return packFloatx80(0, 0, 0);
2079: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low);
2080: }
2081:
2082: if (aExp == 0 && aSig == 0) {
2083: return packFloatx80(0, one_exp, one_sig);
2084: }
2085:
2086: user_rnd_mode = float_rounding_mode;
2087: user_rnd_prec = floatx80_rounding_precision;
2088: float_rounding_mode = float_round_nearest_even;
2089: floatx80_rounding_precision = 80;
2090:
2091: fp0 = a;
2092:
2093: compact = floatx80_make_compact(aExp, aSig);
2094:
2095: if (compact < 0x3FB98000 || compact > 0x400B9B07) { // |X| > 16480 LOG2/LOG10 or |X| < 2^(-70)
2096: if (compact > 0x3FFF8000) { // |X| > 16480
2097: float_rounding_mode = user_rnd_mode;
2098: floatx80_rounding_precision = user_rnd_prec;
2099:
2100: if (aSign) {
2101: return roundAndPackFloatx80(floatx80_rounding_precision, 0, -0x1000, aSig, 0);
2102: } else {
2103: return roundAndPackFloatx80(floatx80_rounding_precision, 0, 0x8000, aSig, 0);
2104: }
2105: } else { // |X| < 2^(-70)
2106: float_rounding_mode = user_rnd_mode;
2107: floatx80_rounding_precision = user_rnd_prec;
2108:
2109: a = floatx80_add(fp0, float32_to_floatx80(0x3F800000)); // 1 + X
2110:
2111: float_raise(float_flag_inexact);
2112:
2113: return a;
2114: }
2115: } else { // 2^(-70) <= |X| <= 16480 LOG 2 / LOG 10
2116: fp1 = fp0; // X
2117: fp1 = floatx80_mul(fp1, float64_to_floatx80(LIT64(0x406A934F0979A371))); // X*64*LOG10/LOG2
2118: n = floatx80_to_int32(fp1); // N=INT(X*64*LOG10/LOG2)
2119: fp1 = int32_to_floatx80(n);
2120:
2121: j = n & 0x3F;
2122: l = n / 64; // NOTE: this is really arithmetic right shift by 6
2123: if (n < 0 && j) { // arithmetic right shift is division and round towards minus infinity
2124: l--;
2125: }
2126: m = l / 2; // NOTE: this is really arithmetic right shift by 1
2127: if (l < 0 && (l & 1)) { // arithmetic right shift is division and round towards minus infinity
2128: m--;
2129: }
2130: m1 = l - m;
2131: m1 += 0x3FFF; // ADJFACT IS 2^(M')
2132:
2133: adjfact = packFloatx80(0, m1, one_sig);
2134: fact1 = exp2_tbl[j];
2135: fact1.high += m;
2136: fact2.high = exp2_tbl2[j]>>16;
2137: fact2.high += m;
2138: fact2.low = (bits64)(exp2_tbl2[j] & 0xFFFF);
2139: fact2.low <<= 48;
2140:
2141: fp2 = fp1; // N
2142: fp1 = floatx80_mul(fp1, float64_to_floatx80(LIT64(0x3F734413509F8000))); // N*(LOG2/64LOG10)_LEAD
2143: fp3 = packFloatx80(1, 0x3FCD, LIT64(0xC0219DC1DA994FD2));
2144: fp2 = floatx80_mul(fp2, fp3); // N*(LOG2/64LOG10)_TRAIL
2145: fp0 = floatx80_sub(fp0, fp1); // X - N L_LEAD
2146: fp0 = floatx80_sub(fp0, fp2); // X - N L_TRAIL
2147: fp2 = packFloatx80(0, 0x4000, LIT64(0x935D8DDDAAA8AC17)); // LOG10
2148: fp0 = floatx80_mul(fp0, fp2); // R
2149:
2150: // EXPR
2151: fp1 = floatx80_mul(fp0, fp0); // S = R*R
2152: fp2 = float64_to_floatx80(LIT64(0x3F56C16D6F7BD0B2)); // A5
2153: fp3 = float64_to_floatx80(LIT64(0x3F811112302C712C)); // A4
2154: fp2 = floatx80_mul(fp2, fp1); // S*A5
2155: fp3 = floatx80_mul(fp3, fp1); // S*A4
2156: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FA5555555554CC1))); // A3+S*A5
2157: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0x3FC5555555554A54))); // A2+S*A4
2158: fp2 = floatx80_mul(fp2, fp1); // S*(A3+S*A5)
2159: fp3 = floatx80_mul(fp3, fp1); // S*(A2+S*A4)
2160: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FE0000000000000))); // A1+S*(A3+S*A5)
2161: fp3 = floatx80_mul(fp3, fp0); // R*S*(A2+S*A4)
2162:
2163: fp2 = floatx80_mul(fp2, fp1); // S*(A1+S*(A3+S*A5))
2164: fp0 = floatx80_add(fp0, fp3); // R+R*S*(A2+S*A4)
2165: fp0 = floatx80_add(fp0, fp2); // EXP(R) - 1
2166:
2167: fp0 = floatx80_mul(fp0, fact1);
2168: fp0 = floatx80_add(fp0, fact2);
2169: fp0 = floatx80_add(fp0, fact1);
2170:
2171: float_rounding_mode = user_rnd_mode;
2172: floatx80_rounding_precision = user_rnd_prec;
2173:
2174: a = floatx80_mul(fp0, adjfact);
2175:
2176: float_raise(float_flag_inexact);
2177:
2178: return a;
2179: }
2180: }
2181:
2182: /*----------------------------------------------------------------------------
2183: | 2 to x
2184: *----------------------------------------------------------------------------*/
2185:
2186: floatx80 floatx80_twotox(floatx80 a)
2187: {
2188: flag aSign;
2189: int32 aExp;
2190: bits64 aSig;
2191:
2192: int8 user_rnd_mode, user_rnd_prec;
2193:
2194: int32 compact, n, j, l, m, m1;
2195: floatx80 fp0, fp1, fp2, fp3, adjfact, fact1, fact2;
2196:
2197: aSig = extractFloatx80Frac(a);
2198: aExp = extractFloatx80Exp(a);
2199: aSign = extractFloatx80Sign(a);
2200:
2201: if (aExp == 0x7FFF) {
2202: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a);
2203: if (aSign) return packFloatx80(0, 0, 0);
2204: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low);
2205: }
2206:
2207: if (aExp == 0 && aSig == 0) {
2208: return packFloatx80(0, one_exp, one_sig);
2209: }
2210:
2211: user_rnd_mode = float_rounding_mode;
2212: user_rnd_prec = floatx80_rounding_precision;
2213: float_rounding_mode = float_round_nearest_even;
2214: floatx80_rounding_precision = 80;
2215:
2216: fp0 = a;
2217:
2218: compact = floatx80_make_compact(aExp, aSig);
2219:
2220: if (compact < 0x3FB98000 || compact > 0x400D80C0) { // |X| > 16480 or |X| < 2^(-70)
2221: if (compact > 0x3FFF8000) { // |X| > 16480
2222: float_rounding_mode = user_rnd_mode;
2223: floatx80_rounding_precision = user_rnd_prec;
2224:
2225: if (aSign) {
2226: return roundAndPackFloatx80(floatx80_rounding_precision, 0, -0x1000, aSig, 0);
2227: } else {
2228: return roundAndPackFloatx80(floatx80_rounding_precision, 0, 0x8000, aSig, 0);
2229: }
2230: } else { // |X| < 2^(-70)
2231: float_rounding_mode = user_rnd_mode;
2232: floatx80_rounding_precision = user_rnd_prec;
2233:
2234: a = floatx80_add(fp0, float32_to_floatx80(0x3F800000)); // 1 + X
2235:
2236: float_raise(float_flag_inexact);
2237:
2238: return a;
2239: }
2240: } else { // 2^(-70) <= |X| <= 16480
2241: fp1 = fp0; // X
2242: fp1 = floatx80_mul(fp1, float32_to_floatx80(0x42800000)); // X * 64
2243: n = floatx80_to_int32(fp1);
2244: fp1 = int32_to_floatx80(n);
2245: j = n & 0x3F;
2246: l = n / 64; // NOTE: this is really arithmetic right shift by 6
2247: if (n < 0 && j) { // arithmetic right shift is division and round towards minus infinity
2248: l--;
2249: }
2250: m = l / 2; // NOTE: this is really arithmetic right shift by 1
2251: if (l < 0 && (l & 1)) { // arithmetic right shift is division and round towards minus infinity
2252: m--;
2253: }
2254: m1 = l - m;
2255: m1 += 0x3FFF; // ADJFACT IS 2^(M')
2256:
2257: adjfact = packFloatx80(0, m1, one_sig);
2258: fact1 = exp2_tbl[j];
2259: fact1.high += m;
2260: fact2.high = exp2_tbl2[j]>>16;
2261: fact2.high += m;
2262: fact2.low = (bits64)(exp2_tbl2[j] & 0xFFFF);
2263: fact2.low <<= 48;
2264:
2265: fp1 = floatx80_mul(fp1, float32_to_floatx80(0x3C800000)); // (1/64)*N
2266: fp0 = floatx80_sub(fp0, fp1); // X - (1/64)*INT(64 X)
2267: fp2 = packFloatx80(0, 0x3FFE, LIT64(0xB17217F7D1CF79AC)); // LOG2
2268: fp0 = floatx80_mul(fp0, fp2); // R
2269:
2270: // EXPR
2271: fp1 = floatx80_mul(fp0, fp0); // S = R*R
2272: fp2 = float64_to_floatx80(LIT64(0x3F56C16D6F7BD0B2)); // A5
2273: fp3 = float64_to_floatx80(LIT64(0x3F811112302C712C)); // A4
2274: fp2 = floatx80_mul(fp2, fp1); // S*A5
2275: fp3 = floatx80_mul(fp3, fp1); // S*A4
2276: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FA5555555554CC1))); // A3+S*A5
2277: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0x3FC5555555554A54))); // A2+S*A4
2278: fp2 = floatx80_mul(fp2, fp1); // S*(A3+S*A5)
2279: fp3 = floatx80_mul(fp3, fp1); // S*(A2+S*A4)
2280: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FE0000000000000))); // A1+S*(A3+S*A5)
2281: fp3 = floatx80_mul(fp3, fp0); // R*S*(A2+S*A4)
2282:
2283: fp2 = floatx80_mul(fp2, fp1); // S*(A1+S*(A3+S*A5))
2284: fp0 = floatx80_add(fp0, fp3); // R+R*S*(A2+S*A4)
2285: fp0 = floatx80_add(fp0, fp2); // EXP(R) - 1
2286:
2287: fp0 = floatx80_mul(fp0, fact1);
2288: fp0 = floatx80_add(fp0, fact2);
2289: fp0 = floatx80_add(fp0, fact1);
2290:
2291: float_rounding_mode = user_rnd_mode;
2292: floatx80_rounding_precision = user_rnd_prec;
2293:
2294: a = floatx80_mul(fp0, adjfact);
2295:
2296: float_raise(float_flag_inexact);
2297:
2298: return a;
2299: }
2300: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.