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