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