![[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:
1.1.1.2 ! root 36: static int32 floatx80_make_compact(int32 aExp, bits64 aSig)
1.1 root 37: {
38: return (aExp<<16)|(aSig>>48);
39: }
40:
41:
42: /*----------------------------------------------------------------------------
43: | Arc cosine
44: *----------------------------------------------------------------------------*/
45:
1.1.1.2 ! root 46: floatx80 floatx80_acos(floatx80 a, float_ctrl* c)
1.1 root 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)) {
1.1.1.2 ! root 62: return propagateFloatx80NaNOneArg(a, c);
1.1 root 63: }
64: if (aExp == 0 && aSig == 0) {
1.1.1.2 ! root 65: float_raise(float_flag_inexact, c);
! 66: return roundAndPackFloatx80(get_float_rounding_precision(c), 0, piby2_exp, pi_sig, 0, c);
1.1 root 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);
1.1.1.2 ! root 75: float_raise(float_flag_inexact, c);
! 76: return floatx80_move(a, c);
1.1 root 77: } else { // X == +1
78: return packFloatx80(0, 0, 0);
79: }
80: } else { // |X| > 1
1.1.1.2 ! root 81: float_raise(float_flag_invalid, c);
1.1 root 82: a.low = floatx80_default_nan_low;
83: a.high = floatx80_default_nan_high;
84: return a;
85: }
86: } // |X| < 1
87:
1.1.1.2 ! root 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);
1.1 root 92:
93: one = packFloatx80(0, one_exp, one_sig);
94: fp0 = a;
95:
1.1.1.2 ! root 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)))
1.1 root 101:
1.1.1.2 ! root 102: set_float_rounding_mode(user_rnd_mode, c);
! 103: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 104:
1.1.1.2 ! root 105: a = floatx80_add(fp0, fp0, c); // 2 * ATAN(SQRT((1-X)/(1+X)))
1.1 root 106:
1.1.1.2 ! root 107: float_raise(float_flag_inexact, c);
1.1 root 108:
109: return a;
110: }
111:
112: /*----------------------------------------------------------------------------
113: | Arc sine
114: *----------------------------------------------------------------------------*/
115:
1.1.1.2 ! root 116: floatx80 floatx80_asin(floatx80 a, float_ctrl* c)
1.1 root 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)) {
1.1.1.2 ! root 132: return propagateFloatx80NaNOneArg(a, c);
1.1 root 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
1.1.1.2 ! root 143: float_raise(float_flag_inexact, c);
1.1 root 144: a = packFloatx80(aSign, piby2_exp, pi_sig);
1.1.1.2 ! root 145: return floatx80_move(a, c);
1.1 root 146: } else { // |X| > 1
1.1.1.2 ! root 147: float_raise(float_flag_invalid, c);
1.1 root 148: a.low = floatx80_default_nan_low;
149: a.high = floatx80_default_nan_high;
150: return a;
151: }
152:
153: } // |X| < 1
154:
1.1.1.2 ! root 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);
1.1 root 159:
160: one = packFloatx80(0, one_exp, one_sig);
161: fp0 = a;
162:
1.1.1.2 ! root 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))
1.1 root 168:
1.1.1.2 ! root 169: set_float_rounding_mode(user_rnd_mode, c);
! 170: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 171:
1.1.1.2 ! root 172: a = floatx80_atan(fp0, c); // ATAN(X/SQRT((1+X)*(1-X)))
1.1 root 173:
1.1.1.2 ! root 174: float_raise(float_flag_inexact, c);
1.1 root 175:
176: return a;
177: }
178:
179: /*----------------------------------------------------------------------------
180: | Arc tangent
181: *----------------------------------------------------------------------------*/
182:
1.1.1.2 ! root 183: floatx80 floatx80_atan(floatx80 a, float_ctrl* c)
1.1 root 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) {
1.1.1.2 ! root 199: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a, c);
1.1 root 200: a = packFloatx80(aSign, piby2_exp, pi_sig);
1.1.1.2 ! root 201: float_raise(float_flag_inexact, c);
! 202: return floatx80_move(a, c);
1.1 root 203: }
204:
205: if (aExp == 0 && aSig == 0) {
206: return packFloatx80(aSign, 0, 0);
207: }
208:
209: compact = floatx80_make_compact(aExp, aSig);
210:
1.1.1.2 ! root 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);
1.1 root 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:
1.1.1.2 ! root 222: set_float_rounding_mode(user_rnd_mode, c);
! 223: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 224:
1.1.1.2 ! root 225: a = floatx80_sub(fp0, fp1, c);
1.1 root 226:
1.1.1.2 ! root 227: float_raise(float_flag_inexact, c);
1.1 root 228:
229: return a;
230: } else {
231: fp0 = a;
232: fp1 = packFloatx80(1, one_exp, one_sig); // -1
1.1.1.2 ! root 233: fp1 = floatx80_div(fp1, fp0, c); // X' = -1/X
1.1 root 234: xsave = fp1;
1.1.1.2 ! root 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);
1.1 root 250: fp1 = packFloatx80(aSign, piby2_exp, pi_sig);
251:
1.1.1.2 ! root 252: set_float_rounding_mode(user_rnd_mode, c);
! 253: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 254:
1.1.1.2 ! root 255: a = floatx80_add(fp0, fp1, c);
1.1 root 256:
1.1.1.2 ! root 257: float_raise(float_flag_inexact, c);
1.1 root 258:
259: return a;
260: }
261: } else { // |X| < 1/16
262: if (compact < 0x3FD78000) { // |X| < 2^(-40)
1.1.1.2 ! root 263: set_float_rounding_mode(user_rnd_mode, c);
! 264: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 265:
1.1.1.2 ! root 266: a = floatx80_move(a, c);
1.1 root 267:
1.1.1.2 ! root 268: float_raise(float_flag_inexact, c);
1.1 root 269:
270: return a;
271: } else {
272: fp0 = a;
273: xsave = a;
1.1.1.2 ! root 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))])
1.1 root 290:
1.1.1.2 ! root 291: set_float_rounding_mode(user_rnd_mode, c);
! 292: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 293:
1.1.1.2 ! root 294: a = floatx80_add(fp0, xsave, c);
1.1 root 295:
1.1.1.2 ! root 296: float_raise(float_flag_inexact, c);
1.1 root 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
1.1.1.2 ! root 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)
1.1 root 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:
1.1.1.2 ! root 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)
1.1 root 334:
1.1.1.2 ! root 335: set_float_rounding_mode(user_rnd_mode, c);
! 336: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 337:
1.1.1.2 ! root 338: a = floatx80_add(fp0, fp3, c); // ATAN(X)
1.1 root 339:
1.1.1.2 ! root 340: float_raise(float_flag_inexact, c);
1.1 root 341:
342: return a;
343: }
344: }
345:
346: /*----------------------------------------------------------------------------
347: | Hyperbolic arc tangent
348: *----------------------------------------------------------------------------*/
349:
1.1.1.2 ! root 350: floatx80 floatx80_atanh(floatx80 a, float_ctrl* c)
1.1 root 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)) {
1.1.1.2 ! root 366: return propagateFloatx80NaNOneArg(a, c);
1.1 root 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
1.1.1.2 ! root 377: float_raise(float_flag_divbyzero, c);
1.1 root 378: return packFloatx80(aSign, 0x7FFF, floatx80_default_infinity_low);
379: } else { // |X| > 1
1.1.1.2 ! root 380: float_raise(float_flag_invalid, c);
1.1 root 381: a.low = floatx80_default_nan_low;
382: a.high = floatx80_default_nan_high;
383: return a;
384: }
385: } // |X| < 1
386:
1.1.1.2 ! root 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);
1.1 root 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
1.1.1.2 ! root 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)
1.1 root 400:
1.1.1.2 ! root 401: set_float_rounding_mode(user_rnd_mode, c);
! 402: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 403:
1.1.1.2 ! root 404: a = floatx80_mul(fp0, fp2, c); // ATANH(X) = SIGN(X) * (1/2) * LOG1P(Z)
1.1 root 405:
1.1.1.2 ! root 406: float_raise(float_flag_inexact, c);
1.1 root 407:
408: return a;
409: }
410:
411: /*----------------------------------------------------------------------------
412: | Cosine
413: *----------------------------------------------------------------------------*/
414:
1.1.1.2 ! root 415: floatx80 floatx80_cos(floatx80 a, float_ctrl* c)
1.1 root 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) {
1.1.1.2 ! root 433: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a, c);
! 434: float_raise(float_flag_invalid, c);
1.1 root 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:
1.1.1.2 ! root 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);
1.1 root 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));
1.1.1.2 ! root 462: fp0 = floatx80_add(fp0, twopi1, c);
1.1 root 463: fp1 = fp0;
1.1.1.2 ! root 464: fp0 = floatx80_add(fp0, twopi2, c);
! 465: fp1 = floatx80_sub(fp1, fp0, c);
! 466: fp1 = floatx80_add(fp1, twopi2, c);
1.1 root 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:
1.1.1.2 ! root 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
1.1 root 495: fp3 = fp0; // FP3 is A
1.1.1.2 ! root 496: fp1 = floatx80_sub(fp1, fp4, c); // FP1 is a := r - p
! 497: fp0 = floatx80_add(fp0, fp1, c); // FP0 is R := A+a
1.1 root 498:
499: if (endflag > 0) {
1.1.1.2 ! root 500: n = floatx80_to_int32(fp2, c);
1.1 root 501: goto sincont;
502: }
1.1.1.2 ! root 503: fp3 = floatx80_sub(fp3, fp0, c); // A-R
! 504: fp1 = floatx80_add(fp1, fp3, c); // FP1 is r := (A-R)+a
1.1 root 505: goto loop;
506: } else {
507: // SINSM
1.1.1.2 ! root 508: fp0 = float32_to_floatx80(0x3F800000, c); // 1
1.1 root 509:
1.1.1.2 ! root 510: set_float_rounding_mode(user_rnd_mode, c);
! 511: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 512:
513: if (adjn) {
514: // COSTINY
1.1.1.2 ! root 515: a = floatx80_sub(fp0, float32_to_floatx80(0x00800000, c), c);
1.1 root 516: } else {
517: // SINTINY
1.1.1.2 ! root 518: a = floatx80_move(a, c);
1.1 root 519: }
1.1.1.2 ! root 520: float_raise(float_flag_inexact, c);
1.1 root 521:
522: return a;
523: }
524: } else {
1.1.1.2 ! root 525: fp1 = floatx80_mul(fp0, float64_to_floatx80(LIT64(0x3FE45F306DC9C883), c), c); // X*2/PI
1.1 root 526:
1.1.1.2 ! root 527: n = floatx80_to_int32(fp1, c);
1.1 root 528: j = 32 + n;
529:
1.1.1.2 ! root 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
1.1 root 532:
533: sincont:
534: if ((n + adjn) & 1) {
535: // COSPOLY
1.1.1.2 ! root 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
1.1 root 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:
1.1.1.2 ! root 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)
1.1 root 560: fp4 = packFloatx80(1, 0x3FF5, LIT64(0xB60B60B60B61D438));
1.1.1.2 ! root 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))
1.1 root 564: fp4 = packFloatx80(0, 0x3FFA, LIT64(0xAAAAAAAAAAAAAB5E));
1.1.1.2 ! root 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)))]
1.1 root 569:
570: x = packFloatx80(xSign, xExp, xSig);
1.1.1.2 ! root 571: fp0 = floatx80_mul(fp0, x, c);
1.1 root 572:
1.1.1.2 ! root 573: set_float_rounding_mode(user_rnd_mode, c);
! 574: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 575:
1.1.1.2 ! root 576: a = floatx80_add(fp0, float32_to_floatx80(posneg1, c), c);
1.1 root 577:
1.1.1.2 ! root 578: float_raise(float_flag_inexact, c);
1.1 root 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:
1.1.1.2 ! root 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)
1.1 root 600: fp4 = packFloatx80(0, 0x3FF8, LIT64(0x88888888888859AF));
1.1.1.2 ! root 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))
1.1 root 604: fp4 = packFloatx80(1, 0x3FFC, LIT64(0xAAAAAAAAAAAAAA99));
1.1.1.2 ! root 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))]
1.1 root 607:
608: x = packFloatx80(xSign, xExp, xSig);
1.1.1.2 ! root 609: fp0 = floatx80_mul(fp0, x, c); // R'*S
! 610: fp0 = floatx80_mul(fp0, fp1, c); // SIN(R')-R'
1.1 root 611:
1.1.1.2 ! root 612: set_float_rounding_mode(user_rnd_mode, c);
! 613: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 614:
1.1.1.2 ! root 615: a = floatx80_add(fp0, x, c);
1.1 root 616:
1.1.1.2 ! root 617: float_raise(float_flag_inexact, c);
1.1 root 618:
619: return a;
620: }
621: }
622: }
623:
624: /*----------------------------------------------------------------------------
625: | Hyperbolic cosine
626: *----------------------------------------------------------------------------*/
627:
1.1.1.2 ! root 628: floatx80 floatx80_cosh(floatx80 a, float_ctrl* c)
1.1 root 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) {
1.1.1.2 ! root 642: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a, c);
1.1 root 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:
1.1.1.2 ! root 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);
1.1 root 654:
655: compact = floatx80_make_compact(aExp, aSig);
656:
657: if (compact > 0x400CB167) {
658: if (compact > 0x400CB2B3) {
1.1.1.2 ! root 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);
1.1 root 662: } else {
663: fp0 = packFloatx80(0, aExp, aSig);
1.1.1.2 ! root 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);
1.1 root 667: fp1 = packFloatx80(0, 0x7FFB, one_sig);
668:
1.1.1.2 ! root 669: set_float_rounding_mode(user_rnd_mode, c);
! 670: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 671:
1.1.1.2 ! root 672: a = floatx80_mul(fp0, fp1, c);
1.1 root 673:
1.1.1.2 ! root 674: float_raise(float_flag_inexact, c);
1.1 root 675:
676: return a;
677: }
678: }
679:
680: fp0 = packFloatx80(0, aExp, aSig); // |X|
1.1.1.2 ! root 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|))
1.1 root 685:
1.1.1.2 ! root 686: set_float_rounding_mode(user_rnd_mode, c);
! 687: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 688:
1.1.1.2 ! root 689: a = floatx80_add(fp0, fp1, c);
1.1 root 690:
1.1.1.2 ! root 691: float_raise(float_flag_inexact, c);
1.1 root 692:
693: return a;
694: }
695:
696: /*----------------------------------------------------------------------------
697: | e to x
698: *----------------------------------------------------------------------------*/
699:
1.1.1.2 ! root 700: floatx80 floatx80_etox(floatx80 a, float_ctrl* c)
1.1 root 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) {
1.1.1.2 ! root 717: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a, c);
1.1 root 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:
1.1.1.2 ! root 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);
1.1 root 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;
1.1.1.2 ! root 739: fp0 = floatx80_mul(fp0, float32_to_floatx80(0x42B8AA3B, c), c); // 64/log2 * X
1.1 root 740: adjflag = 0;
1.1.1.2 ! root 741: n = floatx80_to_int32(fp0, c); // int(64/log2*X)
1.1 root 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
1.1.1.2 ! root 753: fp0 = floatx80_mul(fp0, float32_to_floatx80(0xBC317218, c), c); // N * L1, L1 = lead(-log2/64)
1.1 root 754: l2 = packFloatx80(0, 0x3FDC, LIT64(0x82E308654361C4C6));
1.1.1.2 ! root 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
1.1 root 772:
773: fp1 = exp_tbl[j];
1.1.1.2 ! root 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)
1.1 root 777:
778: scale = packFloatx80(0, m, one_sig);
779: if (adjflag) {
780: adjscale = packFloatx80(0, m1, one_sig);
1.1.1.2 ! root 781: fp0 = floatx80_mul(fp0, adjscale, c);
1.1 root 782: }
783:
1.1.1.2 ! root 784: set_float_rounding_mode(user_rnd_mode, c);
! 785: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 786:
1.1.1.2 ! root 787: a = floatx80_mul(fp0, scale, c);
1.1 root 788:
1.1.1.2 ! root 789: float_raise(float_flag_inexact, c);
1.1 root 790:
791: return a;
792: } else { // |X| >= 16380 log2
793: if (compact > 0x400CB27C) { // |X| >= 16480 log2
1.1.1.2 ! root 794: set_float_rounding_mode(user_rnd_mode, c);
! 795: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 796: if (aSign) {
1.1.1.2 ! root 797: a = roundAndPackFloatx80(get_float_rounding_precision(c), 0, -0x1000, aSig, 0, c);
1.1 root 798: } else {
1.1.1.2 ! root 799: a = roundAndPackFloatx80(get_float_rounding_precision(c), 0, 0x8000, aSig, 0, c);
1.1 root 800: }
1.1.1.2 ! root 801: float_raise(float_flag_inexact, c);
1.1 root 802:
803: return a;
804: } else {
805: fp0 = a;
806: fp1 = a;
1.1.1.2 ! root 807: fp0 = floatx80_mul(fp0, float32_to_floatx80(0x42B8AA3B, c), c); // 64/log2 * X
1.1 root 808: adjflag = 1;
1.1.1.2 ! root 809: n = floatx80_to_int32(fp0, c); // int(64/log2*X)
1.1 root 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)
1.1.1.2 ! root 829: set_float_rounding_mode(user_rnd_mode, c);
! 830: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 831:
1.1.1.2 ! root 832: a = floatx80_add(a, float32_to_floatx80(0x3F800000, c), c); // 1 + X
1.1 root 833:
1.1.1.2 ! root 834: float_raise(float_flag_inexact, c);
1.1 root 835:
836: return a;
837: }
838: }
839:
840: /*----------------------------------------------------------------------------
841: | e to x minus 1
842: *----------------------------------------------------------------------------*/
843:
1.1.1.2 ! root 844: floatx80 floatx80_etoxm1(floatx80 a, float_ctrl* c)
1.1 root 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) {
1.1.1.2 ! root 860: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a, c);
1.1 root 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:
1.1.1.2 ! root 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);
1.1 root 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;
1.1.1.2 ! root 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)
1.1 root 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
1.1.1.2 ! root 894: fp0 = floatx80_mul(fp0, float32_to_floatx80(0xBC317218, c), c); // N * L1, L1 = lead(-log2/64)
1.1 root 895: l2 = packFloatx80(0, 0x3FDC, LIT64(0x82E308654361C4C6));
1.1.1.2 ! root 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
1.1 root 915:
1.1.1.2 ! root 916: fp0 = floatx80_mul(fp0, exp_tbl[j], c); // 2^(J/64)*(Exp(R)-1)
1.1 root 917:
918: if (m >= 64) {
1.1.1.2 ! root 919: fp1 = float32_to_floatx80(exp_tbl2[j], c);
1.1 root 920: onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); // -2^(-M)
1.1.1.2 ! root 921: fp1 = floatx80_add(fp1, onebysc, c);
! 922: fp0 = floatx80_add(fp0, fp1, c);
! 923: fp0 = floatx80_add(fp0, exp_tbl[j], c);
1.1 root 924: } else if (m < -3) {
1.1.1.2 ! root 925: fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j], c), c);
! 926: fp0 = floatx80_add(fp0, exp_tbl[j], c);
1.1 root 927: onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); // -2^(-M)
1.1.1.2 ! root 928: fp0 = floatx80_add(fp0, onebysc, c);
1.1 root 929: } else { // -3 <= m <= 63
930: fp1 = exp_tbl[j];
1.1.1.2 ! root 931: fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j], c), c);
1.1 root 932: onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); // -2^(-M)
1.1.1.2 ! root 933: fp1 = floatx80_add(fp1, onebysc, c);
! 934: fp0 = floatx80_add(fp0, fp1, c);
1.1 root 935: }
936:
937: sc = packFloatx80(0, m + 0x3FFF, one_sig);
938:
1.1.1.2 ! root 939: set_float_rounding_mode(user_rnd_mode, c);
! 940: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 941:
1.1.1.2 ! root 942: a = floatx80_mul(fp0, sc, c);
1.1 root 943:
1.1.1.2 ! root 944: float_raise(float_flag_inexact, c);
1.1 root 945:
946: return a;
947: } else { // |X| > 70 log2
948: if (aSign) {
1.1.1.2 ! root 949: fp0 = float32_to_floatx80(0xBF800000, c); // -1
1.1 root 950:
1.1.1.2 ! root 951: set_float_rounding_mode(user_rnd_mode, c);
! 952: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 953:
1.1.1.2 ! root 954: a = floatx80_add(fp0, float32_to_floatx80(0x00800000, c), c); // -1 + 2^(-126)
1.1 root 955:
1.1.1.2 ! root 956: float_raise(float_flag_inexact, c);
1.1 root 957:
958: return a;
959: } else {
1.1.1.2 ! root 960: set_float_rounding_mode(user_rnd_mode, c);
! 961: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 962:
1.1.1.2 ! root 963: return floatx80_etox(a, c);
1.1 root 964: }
965: }
966: } else { // |X| < 1/4
967: if (aExp >= 0x3FBE) {
968: fp0 = a;
1.1.1.2 ! root 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
1.1 root 989: fp3 = packFloatx80(0, 0x3FFC, LIT64(0xAAAAAAAAAAAAAAAB));
1.1.1.2 ! root 990: fp1 = floatx80_add(fp1, fp3, c); // B2
! 991: fp2 = floatx80_mul(fp2, fp0, c);
! 992: fp1 = floatx80_mul(fp1, fp0, c);
1.1 root 993:
1.1.1.2 ! root 994: fp2 = floatx80_mul(fp2, fp0, c);
! 995: fp1 = floatx80_mul(fp1, a, c);
1.1 root 996:
1.1.1.2 ! root 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
1.1 root 1000:
1.1.1.2 ! root 1001: set_float_rounding_mode(user_rnd_mode, c);
! 1002: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1003:
1.1.1.2 ! root 1004: a = floatx80_add(fp0, a, c);
1.1 root 1005:
1.1.1.2 ! root 1006: float_raise(float_flag_inexact, c);
1.1 root 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)
1.1.1.2 ! root 1014: fp0 = floatx80_mul(fp0, float64_to_floatx80(LIT64(0x48B0000000000000), c), c);
! 1015: fp0 = floatx80_add(fp0, sc, c);
1.1 root 1016:
1.1.1.2 ! root 1017: set_float_rounding_mode(user_rnd_mode, c);
! 1018: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1019:
1.1.1.2 ! root 1020: a = floatx80_mul(fp0, float64_to_floatx80(LIT64(0x3730000000000000), c), c);
1.1 root 1021: } else {
1.1.1.2 ! root 1022: set_float_rounding_mode(user_rnd_mode, c);
! 1023: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1024:
1.1.1.2 ! root 1025: a = floatx80_add(fp0, sc, c);
1.1 root 1026: }
1027:
1.1.1.2 ! root 1028: float_raise(float_flag_inexact, c);
1.1 root 1029:
1030: return a;
1031: }
1032: }
1033: }
1034:
1035: /*----------------------------------------------------------------------------
1036: | Log base 10
1037: *----------------------------------------------------------------------------*/
1038:
1.1.1.2 ! root 1039: floatx80 floatx80_log10(floatx80 a, float_ctrl* c)
1.1 root 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) {
1.1.1.2 ! root 1054: if ((bits64) (aSig<<1)) propagateFloatx80NaNOneArg(a, c);
1.1 root 1055: if (aSign == 0)
1056: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low);
1057: }
1058:
1059: if (aExp == 0 && aSig == 0) {
1.1.1.2 ! root 1060: float_raise(float_flag_divbyzero, c);
1.1 root 1061: return packFloatx80(1, 0x7FFF, floatx80_default_infinity_low);
1062: }
1063:
1064: if (aSign) {
1.1.1.2 ! root 1065: float_raise(float_flag_invalid, c);
1.1 root 1066: a.low = floatx80_default_nan_low;
1067: a.high = floatx80_default_nan_high;
1068: return a;
1069: }
1070:
1.1.1.2 ! root 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);
1.1 root 1075:
1.1.1.2 ! root 1076: fp0 = floatx80_logn(a, c);
1.1 root 1077: fp1 = packFloatx80(0, 0x3FFD, LIT64(0xDE5BD8A937287195)); // INV_L10
1078:
1.1.1.2 ! root 1079: set_float_rounding_mode(user_rnd_mode, c);
! 1080: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1081:
1.1.1.2 ! root 1082: a = floatx80_mul(fp0, fp1, c); // LOGN(X)*INV_L10
1.1 root 1083:
1.1.1.2 ! root 1084: float_raise(float_flag_inexact, c);
1.1 root 1085:
1086: return a;
1087: }
1088:
1089: /*----------------------------------------------------------------------------
1090: | Log base 2
1091: *----------------------------------------------------------------------------*/
1092:
1.1.1.2 ! root 1093: floatx80 floatx80_log2(floatx80 a, float_ctrl* c)
1.1 root 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) {
1.1.1.2 ! root 1108: if ((bits64) (aSig<<1)) propagateFloatx80NaNOneArg(a, c);
1.1 root 1109: if (aSign == 0)
1110: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low);
1111: }
1112:
1113: if (aExp == 0) {
1114: if (aSig == 0) {
1.1.1.2 ! root 1115: float_raise(float_flag_divbyzero, c);
1.1 root 1116: return packFloatx80(1, 0x7FFF, floatx80_default_infinity_low);
1117: }
1118: normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
1119: }
1120:
1121: if (aSign) {
1.1.1.2 ! root 1122: float_raise(float_flag_invalid, c);
1.1 root 1123: a.low = floatx80_default_nan_low;
1124: a.high = floatx80_default_nan_high;
1125: return a;
1126: }
1127:
1.1.1.2 ! root 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);
1.1 root 1132:
1133: if (aSig == one_sig) { // X is 2^k
1.1.1.2 ! root 1134: set_float_rounding_mode(user_rnd_mode, c);
! 1135: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1136:
1137: a = int32_to_floatx80(aExp-0x3FFF);
1138: } else {
1.1.1.2 ! root 1139: fp0 = floatx80_logn(a, c);
1.1 root 1140: fp1 = packFloatx80(0, 0x3FFF, LIT64(0xB8AA3B295C17F0BC)); // INV_L2
1141:
1.1.1.2 ! root 1142: set_float_rounding_mode(user_rnd_mode, c);
! 1143: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1144:
1.1.1.2 ! root 1145: a = floatx80_mul(fp0, fp1, c); // LOGN(X)*INV_L2
1.1 root 1146: }
1147:
1.1.1.2 ! root 1148: float_raise(float_flag_inexact, c);
1.1 root 1149:
1150: return a;
1151: }
1152:
1153: /*----------------------------------------------------------------------------
1154: | Log base e
1155: *----------------------------------------------------------------------------*/
1156:
1.1.1.2 ! root 1157: floatx80 floatx80_logn(floatx80 a, float_ctrl* c)
1.1 root 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) {
1.1.1.2 ! root 1173: if ((bits64) (aSig<<1)) propagateFloatx80NaNOneArg(a, c);
1.1 root 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
1.1.1.2 ! root 1182: float_raise(float_flag_divbyzero, c);
1.1 root 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) {
1.1.1.2 ! root 1198: float_raise(float_flag_invalid, c);
1.1 root 1199: a.low = floatx80_default_nan_low;
1200: a.high = floatx80_default_nan_high;
1201: return a;
1202: }
1203:
1.1.1.2 ! root 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);
1.1 root 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:
1.1.1.2 ! root 1222: fp0 = floatx80_sub(fp0, f, c); // Y-F
1.1 root 1223:
1224: // LP1CONT1
1.1.1.2 ! root 1225: fp0 = floatx80_mul(fp0, log_tbl[j], c); // FP0 IS U = (Y-F)/F
1.1 root 1226: logof2 = packFloatx80(0, 0x3FFE, LIT64(0xB17217F7D1CF79AC));
1.1.1.2 ! root 1227: klog2 = floatx80_mul(fp1, logof2, c); // FP1 IS K*LOG2
! 1228: fp2 = floatx80_mul(fp0, fp0, c); // FP2 IS V=U*U
1.1 root 1229:
1230: fp3 = fp2;
1231: fp1 = fp2;
1232:
1.1.1.2 ! root 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)
1.1 root 1248:
1.1.1.2 ! root 1249: set_float_rounding_mode(user_rnd_mode, c);
! 1250: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1251:
1.1.1.2 ! root 1252: a = floatx80_add(fp0, klog2, c);
1.1 root 1253:
1.1.1.2 ! root 1254: float_raise(float_flag_inexact, c);
1.1 root 1255:
1256: return a;
1257: } else { // |X-1| >= 1/16
1258: fp0 = a;
1259: fp1 = a;
1.1.1.2 ! root 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)
1.1 root 1263:
1264: // LP1CONT2
1.1.1.2 ! root 1265: fp1 = floatx80_div(fp1, fp0, c); // U
1.1 root 1266: saveu = fp1;
1.1.1.2 ! root 1267: fp0 = floatx80_mul(fp1, fp1, c); // FP0 IS V = U*U
! 1268: fp1 = floatx80_mul(fp0, fp0, c); // FP1 IS W = V*V
1.1 root 1269:
1.1.1.2 ! root 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)] )
1.1 root 1283:
1.1.1.2 ! root 1284: set_float_rounding_mode(user_rnd_mode, c);
! 1285: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1286:
1.1.1.2 ! root 1287: a = floatx80_add(fp0, saveu, c);
1.1 root 1288:
1289: //if (!floatx80_is_zero(a)) {
1.1.1.2 ! root 1290: float_raise(float_flag_inexact, c);
1.1 root 1291: //}
1292:
1293: return a;
1294: }
1295: }
1296:
1297: /*----------------------------------------------------------------------------
1298: | Log base e of x plus 1
1299: *----------------------------------------------------------------------------*/
1300:
1.1.1.2 ! root 1301: floatx80 floatx80_lognp1(floatx80 a, float_ctrl* c)
1.1 root 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) {
1.1.1.2 ! root 1317: if ((bits64) (aSig<<1)) propagateFloatx80NaNOneArg(a, c);
1.1 root 1318: if (aSign) {
1.1.1.2 ! root 1319: float_raise(float_flag_invalid, c);
1.1 root 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) {
1.1.1.2 ! root 1333: float_raise(float_flag_divbyzero, c);
1.1 root 1334: packFloatx80(aSign, 0x7FFF, floatx80_default_infinity_low);
1335: }
1.1.1.2 ! root 1336: float_raise(float_flag_invalid, c);
1.1 root 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
1.1.1.2 ! root 1343: float_raise(float_flag_inexact, c);
! 1344: return floatx80_move(a, c);
1.1 root 1345: }
1346:
1.1.1.2 ! root 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);
1.1 root 1351:
1352: compact = floatx80_make_compact(aExp, aSig);
1353:
1354: fp0 = a; // Z
1355: fp1 = a;
1356:
1.1.1.2 ! root 1357: fp0 = floatx80_add(fp0, float32_to_floatx80(0x3F800000, c), c); // X = (1+Z)
1.1 root 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:
1.1.1.2 ! root 1374: fp0 = floatx80_sub(fp0, f, c); // Y-F
1.1 root 1375:
1376: lp1cont1:
1377: // LP1CONT1
1.1.1.2 ! root 1378: fp0 = floatx80_mul(fp0, log_tbl[j], c); // FP0 IS U = (Y-F)/F
1.1 root 1379: logof2 = packFloatx80(0, 0x3FFE, LIT64(0xB17217F7D1CF79AC));
1.1.1.2 ! root 1380: klog2 = floatx80_mul(fp1, logof2, c); // FP1 IS K*LOG2
! 1381: fp2 = floatx80_mul(fp0, fp0, c); // FP2 IS V=U*U
1.1 root 1382:
1383: fp3 = fp2;
1384: fp1 = fp2;
1385:
1.1.1.2 ! root 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)
1.1 root 1401:
1.1.1.2 ! root 1402: set_float_rounding_mode(user_rnd_mode, c);
! 1403: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1404:
1.1.1.2 ! root 1405: a = floatx80_add(fp0, klog2, c);
1.1 root 1406:
1.1.1.2 ! root 1407: float_raise(float_flag_inexact, c);
1.1 root 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
1.1.1.2 ! root 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
1.1 root 1420: fp1 = packFloatx80(0, 0, 0); // K = 0
1421: } else {
1422: // KISNEG
1.1.1.2 ! root 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
1.1 root 1426: fp1 = packFloatx80(1, one_exp, one_sig); // K = -1
1427: }
1428: goto lp1cont1;
1429: } else {
1430: // LP1ONE16
1.1.1.2 ! root 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
1.1 root 1433:
1434: // LP1CONT2
1.1.1.2 ! root 1435: fp1 = floatx80_div(fp1, fp0, c); // U
1.1 root 1436: saveu = fp1;
1.1.1.2 ! root 1437: fp0 = floatx80_mul(fp1, fp1, c); // FP0 IS V = U*U
! 1438: fp1 = floatx80_mul(fp0, fp0, c); // FP1 IS W = V*V
1.1 root 1439:
1.1.1.2 ! root 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)] )
1.1 root 1453:
1.1.1.2 ! root 1454: set_float_rounding_mode(user_rnd_mode, c);
! 1455: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1456:
1.1.1.2 ! root 1457: a = floatx80_add(fp0, saveu, c);
1.1 root 1458:
1459: //if (!floatx80_is_zero(a)) {
1.1.1.2 ! root 1460: float_raise(float_flag_inexact, c);
1.1 root 1461: //}
1462:
1463: return a;
1464: }
1465: }
1466:
1467: /*----------------------------------------------------------------------------
1468: | Sine
1469: *----------------------------------------------------------------------------*/
1470:
1.1.1.2 ! root 1471: floatx80 floatx80_sin(floatx80 a, float_ctrl* c)
1.1 root 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) {
1.1.1.2 ! root 1489: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a, c);
! 1490: float_raise(float_flag_invalid, c);
1.1 root 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:
1.1.1.2 ! root 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);
1.1 root 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));
1.1.1.2 ! root 1518: fp0 = floatx80_add(fp0, twopi1, c);
1.1 root 1519: fp1 = fp0;
1.1.1.2 ! root 1520: fp0 = floatx80_add(fp0, twopi2, c);
! 1521: fp1 = floatx80_sub(fp1, fp0, c);
! 1522: fp1 = floatx80_add(fp1, twopi2, c);
1.1 root 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:
1.1.1.2 ! root 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
1.1 root 1551: fp3 = fp0; // FP3 is A
1.1.1.2 ! root 1552: fp1 = floatx80_sub(fp1, fp4, c); // FP1 is a := r - p
! 1553: fp0 = floatx80_add(fp0, fp1, c); // FP0 is R := A+a
1.1 root 1554:
1555: if (endflag > 0) {
1.1.1.2 ! root 1556: n = floatx80_to_int32(fp2, c);
1.1 root 1557: goto sincont;
1558: }
1.1.1.2 ! root 1559: fp3 = floatx80_sub(fp3, fp0, c); // A-R
! 1560: fp1 = floatx80_add(fp1, fp3, c); // FP1 is r := (A-R)+a
1.1 root 1561: goto loop;
1562: } else {
1563: // SINSM
1.1.1.2 ! root 1564: fp0 = float32_to_floatx80(0x3F800000, c); // 1
1.1 root 1565:
1.1.1.2 ! root 1566: set_float_rounding_mode(user_rnd_mode, c);
! 1567: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1568:
1569: if (adjn) {
1570: // COSTINY
1.1.1.2 ! root 1571: a = floatx80_sub(fp0, float32_to_floatx80(0x00800000, c), c);
1.1 root 1572: } else {
1573: // SINTINY
1.1.1.2 ! root 1574: a = floatx80_move(a, c);
1.1 root 1575: }
1.1.1.2 ! root 1576: float_raise(float_flag_inexact, c);
1.1 root 1577:
1578: return a;
1579: }
1580: } else {
1.1.1.2 ! root 1581: fp1 = floatx80_mul(fp0, float64_to_floatx80(LIT64(0x3FE45F306DC9C883), c), c); // X*2/PI
1.1 root 1582:
1.1.1.2 ! root 1583: n = floatx80_to_int32(fp1, c);
1.1 root 1584: j = 32 + n;
1585:
1.1.1.2 ! root 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
1.1 root 1588:
1589: sincont:
1590: if ((n + adjn) & 1) {
1591: // COSPOLY
1.1.1.2 ! root 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
1.1 root 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:
1.1.1.2 ! root 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)
1.1 root 1616: fp4 = packFloatx80(1, 0x3FF5, LIT64(0xB60B60B60B61D438));
1.1.1.2 ! root 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))
1.1 root 1620: fp4 = packFloatx80(0, 0x3FFA, LIT64(0xAAAAAAAAAAAAAB5E));
1.1.1.2 ! root 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)))]
1.1 root 1625:
1626: x = packFloatx80(xSign, xExp, xSig);
1.1.1.2 ! root 1627: fp0 = floatx80_mul(fp0, x, c);
1.1 root 1628:
1.1.1.2 ! root 1629: set_float_rounding_mode(user_rnd_mode, c);
! 1630: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1631:
1.1.1.2 ! root 1632: a = floatx80_add(fp0, float32_to_floatx80(posneg1, c), c);
1.1 root 1633:
1.1.1.2 ! root 1634: float_raise(float_flag_inexact, c);
1.1 root 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:
1.1.1.2 ! root 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)
1.1 root 1656: fp4 = packFloatx80(0, 0x3FF8, LIT64(0x88888888888859AF));
1.1.1.2 ! root 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))
1.1 root 1660: fp4 = packFloatx80(1, 0x3FFC, LIT64(0xAAAAAAAAAAAAAA99));
1.1.1.2 ! root 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))]
1.1 root 1663:
1664: x = packFloatx80(xSign, xExp, xSig);
1.1.1.2 ! root 1665: fp0 = floatx80_mul(fp0, x, c); // R'*S
! 1666: fp0 = floatx80_mul(fp0, fp1, c); // SIN(R')-R'
1.1 root 1667:
1.1.1.2 ! root 1668: set_float_rounding_mode(user_rnd_mode, c);
! 1669: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1670:
1.1.1.2 ! root 1671: a = floatx80_add(fp0, x, c);
1.1 root 1672:
1.1.1.2 ! root 1673: float_raise(float_flag_inexact, c);
1.1 root 1674:
1675: return a;
1676: }
1677: }
1678: }
1679:
1680: /*----------------------------------------------------------------------------
1681: | Hyperbolic sine
1682: *----------------------------------------------------------------------------*/
1683:
1.1.1.2 ! root 1684: floatx80 floatx80_sinh(floatx80 a, float_ctrl* c)
1.1 root 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) {
1.1.1.2 ! root 1701: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a, c);
1.1 root 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:
1.1.1.2 ! root 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);
1.1 root 1713:
1714: compact = floatx80_make_compact(aExp, aSig);
1715:
1716: if (compact > 0x400CB167) {
1717: // SINHBIG
1718: if (compact > 0x400CB2B3) {
1.1.1.2 ! root 1719: set_float_rounding_mode(user_rnd_mode, c);
! 1720: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1721:
1.1.1.2 ! root 1722: return roundAndPackFloatx80(get_float_rounding_precision(c), aSign, 0x8000, aSig, 0, c);
1.1 root 1723: } else {
1.1.1.2 ! root 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);
1.1 root 1728: fp2 = packFloatx80(aSign, 0x7FFB, one_sig);
1729:
1.1.1.2 ! root 1730: set_float_rounding_mode(user_rnd_mode, c);
! 1731: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1732:
1.1.1.2 ! root 1733: a = floatx80_mul(fp0, fp2, c);
1.1 root 1734:
1.1.1.2 ! root 1735: float_raise(float_flag_inexact, c);
1.1 root 1736:
1737: return a;
1738: }
1739: } else { // |X| < 16380 LOG2
1.1.1.2 ! root 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
1.1 root 1743: fp2 = fp0;
1.1.1.2 ! root 1744: fp0 = floatx80_div(fp0, fp1, c); // Z/(1+Z)
! 1745: fp0 = floatx80_add(fp0, fp2, c);
1.1 root 1746:
1747: fact = 0x3F000000;
1748: fact |= aSign ? 0x80000000 : 0x00000000;
1749:
1.1.1.2 ! root 1750: set_float_rounding_mode(user_rnd_mode, c);
! 1751: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1752:
1.1.1.2 ! root 1753: a = floatx80_mul(fp0, float32_to_floatx80(fact, c), c);
1.1 root 1754:
1.1.1.2 ! root 1755: float_raise(float_flag_inexact, c);
1.1 root 1756:
1757: return a;
1758: }
1759: }
1760:
1761: /*----------------------------------------------------------------------------
1762: | Tangent
1763: *----------------------------------------------------------------------------*/
1764:
1.1.1.2 ! root 1765: floatx80 floatx80_tan(floatx80 a, float_ctrl* c)
1.1 root 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) {
1.1.1.2 ! root 1783: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a, c);
! 1784: float_raise(float_flag_invalid, c);
1.1 root 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:
1.1.1.2 ! root 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);
1.1 root 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));
1.1.1.2 ! root 1810: fp0 = floatx80_add(fp0, twopi1, c);
1.1 root 1811: fp1 = fp0;
1.1.1.2 ! root 1812: fp0 = floatx80_add(fp0, twopi2, c);
! 1813: fp1 = floatx80_sub(fp1, fp0, c);
! 1814: fp1 = floatx80_add(fp1, twopi2, c);
1.1 root 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:
1.1.1.2 ! root 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
1.1 root 1843: fp3 = fp0; // FP3 is A
1.1.1.2 ! root 1844: fp1 = floatx80_sub(fp1, fp4, c); // FP1 is a := r - p
! 1845: fp0 = floatx80_add(fp0, fp1, c); // FP0 is R := A+a
1.1 root 1846:
1847: if (endflag > 0) {
1.1.1.2 ! root 1848: n = floatx80_to_int32(fp2, c);
1.1 root 1849: goto tancont;
1850: }
1.1.1.2 ! root 1851: fp3 = floatx80_sub(fp3, fp0, c); // A-R
! 1852: fp1 = floatx80_add(fp1, fp3, c); // FP1 is r := (A-R)+a
1.1 root 1853: goto loop;
1854: } else {
1.1.1.2 ! root 1855: set_float_rounding_mode(user_rnd_mode, c);
! 1856: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1857:
1.1.1.2 ! root 1858: a = floatx80_move(a, c);
1.1 root 1859:
1.1.1.2 ! root 1860: float_raise(float_flag_inexact, c);
1.1 root 1861:
1862: return a;
1863: }
1864: } else {
1.1.1.2 ! root 1865: fp1 = floatx80_mul(fp0, float64_to_floatx80(LIT64(0x3FE45F306DC9C883), c), c); // X*2/PI
1.1 root 1866:
1.1.1.2 ! root 1867: n = floatx80_to_int32(fp1, c);
1.1 root 1868: j = 32 + n;
1869:
1.1.1.2 ! root 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
1.1 root 1872:
1873: tancont:
1874: if (n & 1) {
1875: // NODD
1876: fp1 = fp0; // R
1.1.1.2 ! root 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
1.1 root 1883: fp4 = packFloatx80(0, 0x3FF6, LIT64(0xE073D3FC199C4A00));
1.1.1.2 ! root 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)
1.1 root 1887: fp4 = packFloatx80(0, 0x3FF9, LIT64(0xD23CD68415D95FA1));
1.1.1.2 ! root 1888: fp3 = floatx80_add(fp3, fp4, c); // Q2+S(Q3+SQ4)
1.1 root 1889: fp4 = packFloatx80(1, 0x3FFC, LIT64(0x8895A6C5FB423BCA));
1.1.1.2 ! root 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))
1.1 root 1893: fp4 = packFloatx80(1, 0x3FFD, LIT64(0xEEF57E0DA84BC8CE));
1.1.1.2 ! root 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)))
1.1 root 1899:
1900: xSign = extractFloatx80Sign(fp1);
1901: xExp = extractFloatx80Exp(fp1);
1902: xSig = extractFloatx80Frac(fp1);
1903: xSign ^= 1;
1904: fp1 = packFloatx80(xSign, xExp, xSig);
1905:
1.1.1.2 ! root 1906: set_float_rounding_mode(user_rnd_mode, c);
! 1907: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1908:
1.1.1.2 ! root 1909: a = floatx80_div(fp0, fp1, c);
1.1 root 1910:
1.1.1.2 ! root 1911: float_raise(float_flag_inexact, c);
1.1 root 1912:
1913: return a;
1914: } else {
1.1.1.2 ! root 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
1.1 root 1921: fp4 = packFloatx80(0, 0x3FF6, LIT64(0xE073D3FC199C4A00));
1.1.1.2 ! root 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)
1.1 root 1925: fp4 = packFloatx80(0, 0x3FF9, LIT64(0xD23CD68415D95FA1));
1.1.1.2 ! root 1926: fp3 = floatx80_add(fp3, fp4, c); // Q2+S(Q3+SQ4)
1.1 root 1927: fp4 = packFloatx80(1, 0x3FFC, LIT64(0x8895A6C5FB423BCA));
1.1.1.2 ! root 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))
1.1 root 1931: fp4 = packFloatx80(1, 0x3FFD, LIT64(0xEEF57E0DA84BC8CE));
1.1.1.2 ! root 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)))
1.1 root 1937:
1.1.1.2 ! root 1938: set_float_rounding_mode(user_rnd_mode, c);
! 1939: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1940:
1.1.1.2 ! root 1941: a = floatx80_div(fp0, fp1, c);
1.1 root 1942:
1.1.1.2 ! root 1943: float_raise(float_flag_inexact, c);
1.1 root 1944:
1945: return a;
1946: }
1947: }
1948: }
1949:
1950: /*----------------------------------------------------------------------------
1951: | Hyperbolic tangent
1952: *----------------------------------------------------------------------------*/
1953:
1.1.1.2 ! root 1954: floatx80 floatx80_tanh(floatx80 a, float_ctrl* c)
1.1 root 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) {
1.1.1.2 ! root 1971: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a, c);
1.1 root 1972: return packFloatx80(aSign, one_exp, one_sig);
1973: }
1974:
1975: if (aExp == 0 && aSig == 0) {
1976: return packFloatx80(aSign, 0, 0);
1977: }
1978:
1.1.1.2 ! root 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);
1.1 root 1983:
1984: compact = floatx80_make_compact(aExp, aSig);
1985:
1986: if (compact < 0x3FD78000 || compact > 0x3FFFDDCE) {
1987: // TANHBORS
1988: if (compact < 0x3FFF8000) {
1989: // TANHSM
1.1.1.2 ! root 1990: set_float_rounding_mode(user_rnd_mode, c);
! 1991: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 1992:
1.1.1.2 ! root 1993: a = floatx80_move(a, c);
1.1 root 1994:
1.1.1.2 ! root 1995: float_raise(float_flag_inexact, c);
1.1 root 1996:
1997: return a;
1998: } else {
1999: if (compact > 0x40048AA1) {
2000: // TANHHUGE
2001: sign = 0x3F800000;
2002: sign |= aSign ? 0x80000000 : 0x00000000;
1.1.1.2 ! root 2003: fp0 = float32_to_floatx80(sign, c);
1.1 root 2004: sign &= 0x80000000;
2005: sign ^= 0x80800000; // -SIGN(X)*EPS
2006:
1.1.1.2 ! root 2007: set_float_rounding_mode(user_rnd_mode, c);
! 2008: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 2009:
1.1.1.2 ! root 2010: a = floatx80_add(fp0, float32_to_floatx80(sign, c), c);
1.1 root 2011:
1.1.1.2 ! root 2012: float_raise(float_flag_inexact, c);
1.1 root 2013:
2014: return a;
2015: } else {
2016: fp0 = packFloatx80(0, aExp+1, aSig); // Y = 2|X|
1.1.1.2 ! root 2017: fp0 = floatx80_etox(fp0, c); // FP0 IS EXP(Y)
! 2018: fp0 = floatx80_add(fp0, float32_to_floatx80(0x3F800000, c), c); // EXP(Y)+1
1.1 root 2019: sign = aSign ? 0x80000000 : 0x00000000;
1.1.1.2 ! root 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
1.1 root 2022:
1.1.1.2 ! root 2023: set_float_rounding_mode(user_rnd_mode, c);
! 2024: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 2025:
1.1.1.2 ! root 2026: a = floatx80_add(fp1, fp0, c);
1.1 root 2027:
1.1.1.2 ! root 2028: float_raise(float_flag_inexact, c);
1.1 root 2029:
2030: return a;
2031: }
2032: }
2033: } else { // 2**(-40) < |X| < (5/2)LOG2
2034: fp0 = packFloatx80(0, aExp+1, aSig); // Y = 2|X|
1.1.1.2 ! root 2035: fp0 = floatx80_etoxm1(fp0, c); // FP0 IS Z = EXPM1(Y)
! 2036: fp1 = floatx80_add(fp0, float32_to_floatx80(0x40000000, c), c); // Z+2
1.1 root 2037:
2038: vSign = extractFloatx80Sign(fp1);
2039: vExp = extractFloatx80Exp(fp1);
2040: vSig = extractFloatx80Frac(fp1);
2041:
2042: fp1 = packFloatx80(vSign ^ aSign, vExp, vSig);
2043:
1.1.1.2 ! root 2044: set_float_rounding_mode(user_rnd_mode, c);
! 2045: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 2046:
1.1.1.2 ! root 2047: a = floatx80_div(fp0, fp1, c);
1.1 root 2048:
1.1.1.2 ! root 2049: float_raise(float_flag_inexact, c);
1.1 root 2050:
2051: return a;
2052: }
2053: }
2054:
2055: /*----------------------------------------------------------------------------
2056: | 10 to x
2057: *----------------------------------------------------------------------------*/
2058:
1.1.1.2 ! root 2059: floatx80 floatx80_tentox(floatx80 a, float_ctrl* c)
1.1 root 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) {
1.1.1.2 ! root 2075: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a, c);
1.1 root 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:
1.1.1.2 ! root 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);
1.1 root 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
1.1.1.2 ! root 2095: set_float_rounding_mode(user_rnd_mode, c);
! 2096: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 2097:
2098: if (aSign) {
1.1.1.2 ! root 2099: return roundAndPackFloatx80(get_float_rounding_precision(c), 0, -0x1000, aSig, 0, c);
1.1 root 2100: } else {
1.1.1.2 ! root 2101: return roundAndPackFloatx80(get_float_rounding_precision(c), 0, 0x8000, aSig, 0, c);
1.1 root 2102: }
2103: } else { // |X| < 2^(-70)
1.1.1.2 ! root 2104: set_float_rounding_mode(user_rnd_mode, c);
! 2105: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 2106:
1.1.1.2 ! root 2107: a = floatx80_add(fp0, float32_to_floatx80(0x3F800000, c), c); // 1 + X
1.1 root 2108:
1.1.1.2 ! root 2109: float_raise(float_flag_inexact, c);
1.1 root 2110:
2111: return a;
2112: }
2113: } else { // 2^(-70) <= |X| <= 16480 LOG 2 / LOG 10
2114: fp1 = fp0; // X
1.1.1.2 ! root 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)
1.1 root 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
1.1.1.2 ! root 2140: fp1 = floatx80_mul(fp1, float64_to_floatx80(LIT64(0x3F734413509F8000), c), c); // N*(LOG2/64LOG10)_LEAD
1.1 root 2141: fp3 = packFloatx80(1, 0x3FCD, LIT64(0xC0219DC1DA994FD2));
1.1.1.2 ! root 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
1.1 root 2145: fp2 = packFloatx80(0, 0x4000, LIT64(0x935D8DDDAAA8AC17)); // LOG10
1.1.1.2 ! root 2146: fp0 = floatx80_mul(fp0, fp2, c); // R
1.1 root 2147:
2148: // EXPR
1.1.1.2 ! root 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);
1.1 root 2168:
1.1.1.2 ! root 2169: set_float_rounding_mode(user_rnd_mode, c);
! 2170: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 2171:
1.1.1.2 ! root 2172: a = floatx80_mul(fp0, adjfact, c);
1.1 root 2173:
1.1.1.2 ! root 2174: float_raise(float_flag_inexact, c);
1.1 root 2175:
2176: return a;
2177: }
2178: }
2179:
2180: /*----------------------------------------------------------------------------
2181: | 2 to x
2182: *----------------------------------------------------------------------------*/
2183:
1.1.1.2 ! root 2184: floatx80 floatx80_twotox(floatx80 a, float_ctrl* c)
1.1 root 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) {
1.1.1.2 ! root 2200: if ((bits64) (aSig<<1)) return propagateFloatx80NaNOneArg(a, c);
1.1 root 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:
1.1.1.2 ! root 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);
1.1 root 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
1.1.1.2 ! root 2220: set_float_rounding_mode(user_rnd_mode, c);
! 2221: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 2222:
2223: if (aSign) {
1.1.1.2 ! root 2224: return roundAndPackFloatx80(get_float_rounding_precision(c), 0, -0x1000, aSig, 0, c);
1.1 root 2225: } else {
1.1.1.2 ! root 2226: return roundAndPackFloatx80(get_float_rounding_precision(c), 0, 0x8000, aSig, 0, c);
1.1 root 2227: }
2228: } else { // |X| < 2^(-70)
1.1.1.2 ! root 2229: set_float_rounding_mode(user_rnd_mode, c);
! 2230: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 2231:
1.1.1.2 ! root 2232: a = floatx80_add(fp0, float32_to_floatx80(0x3F800000, c), c); // 1 + X
1.1 root 2233:
1.1.1.2 ! root 2234: float_raise(float_flag_inexact, c);
1.1 root 2235:
2236: return a;
2237: }
2238: } else { // 2^(-70) <= |X| <= 16480
2239: fp1 = fp0; // X
1.1.1.2 ! root 2240: fp1 = floatx80_mul(fp1, float32_to_floatx80(0x42800000, c), c); // X * 64
! 2241: n = floatx80_to_int32(fp1, c);
1.1 root 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:
1.1.1.2 ! root 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)
1.1 root 2265: fp2 = packFloatx80(0, 0x3FFE, LIT64(0xB17217F7D1CF79AC)); // LOG2
1.1.1.2 ! root 2266: fp0 = floatx80_mul(fp0, fp2, c); // R
1.1 root 2267:
2268: // EXPR
1.1.1.2 ! root 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);
1.1 root 2288:
1.1.1.2 ! root 2289: set_float_rounding_mode(user_rnd_mode, c);
! 2290: set_float_rounding_precision(user_rnd_prec, c);
1.1 root 2291:
1.1.1.2 ! root 2292: a = floatx80_mul(fp0, adjfact, c);
1.1 root 2293:
1.1.1.2 ! root 2294: float_raise(float_flag_inexact, c);
1.1 root 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.