@samitouri / QOSamiQemu
Code Commits Issues Pulls Sessions CI/CD
master
master
c 2,822 lines 104 KB
Raw
1 /*
2 * Ported from a work by Andreas Grabher for Previous, NeXT Computer Emulator,
3 * derived from NetBSD M68040 FPSP functions,
4 * derived from release 2a of the SoftFloat IEC/IEEE Floating-point Arithmetic
5 * Package. Those parts of the code (and some later contributions) are
6 * provided under that license, as detailed below.
7 * It has subsequently been modified by contributors to the QEMU Project,
8 * so some portions are provided under:
9 * the SoftFloat-2a license
10 * the BSD license
11 * GPL-v2-or-later
12 *
13 * Any future contributions to this file will be taken to be licensed under
14 * the Softfloat-2a license unless specifically indicated otherwise.
15 */
16
17 /*
18 * Portions of this work are licensed under the terms of the GNU GPL,
19 * version 2 or later. See the COPYING file in the top-level directory.
20 */
21
22 #include "qemu/osdep.h"
23 #include "softfloat.h"
24 #include "fpu/softfloat-macros.h"
25 #include "softfloat_fpsp_tables.h"
26
27 #define pi_exp 0x4000
28 #define piby2_exp 0x3FFF
29 #define pi_sig UINT64_C(0xc90fdaa22168c235)
30
31 static floatx80 propagateFloatx80NaNOneArg(floatx80 a, float_status *status)
32 {
33 if (floatx80_is_signaling_nan(a, status)) {
34 float_raise(float_flag_invalid, status);
35 a = floatx80_silence_nan(a, status);
36 }
37
38 if (get_default_nan_mode(status)) {
39 return floatx80_default_nan(status);
40 }
41
42 return a;
43 }
44
45 /*
46 * Returns the mantissa of the extended double-precision floating-point
47 * value `a'.
48 */
49
50 floatx80 floatx80_getman(floatx80 a, float_status *status)
51 {
52 bool aSign;
53 int32_t aExp;
54 uint64_t aSig;
55
56 aSig = extractFloatx80Frac(a);
57 aExp = extractFloatx80Exp(a);
58 aSign = extractFloatx80Sign(a);
59
60 if (aExp == 0x7FFF) {
61 if ((uint64_t) (aSig << 1)) {
62 return propagateFloatx80NaNOneArg(a , status);
63 }
64 float_raise(float_flag_invalid , status);
65 return floatx80_default_nan(status);
66 }
67
68 if (aExp == 0) {
69 if (aSig == 0) {
70 return packFloatx80(aSign, 0, 0);
71 }
72 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
73 }
74
75 return roundAndPackFloatx80(get_floatx80_rounding_precision(status),
76 aSign,
77 0x3FFF, aSig, 0, status);
78 }
79
80 /*
81 * Returns the exponent of the extended double-precision floating-point
82 * value `a' as an extended double-precision value.
83 */
84
85 floatx80 floatx80_getexp(floatx80 a, float_status *status)
86 {
87 bool aSign;
88 int32_t aExp;
89 uint64_t aSig;
90
91 aSig = extractFloatx80Frac(a);
92 aExp = extractFloatx80Exp(a);
93 aSign = extractFloatx80Sign(a);
94
95 if (aExp == 0x7FFF) {
96 if ((uint64_t) (aSig << 1)) {
97 return propagateFloatx80NaNOneArg(a , status);
98 }
99 float_raise(float_flag_invalid , status);
100 return floatx80_default_nan(status);
101 }
102
103 if (aExp == 0) {
104 if (aSig == 0) {
105 return packFloatx80(aSign, 0, 0);
106 }
107 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
108 }
109
110 return int32_to_floatx80(aExp - 0x3FFF, status);
111 }
112
113 /*
114 * Scales extended double-precision floating-point value in operand `a' by
115 * value `b'. The function truncates the value in the second operand 'b' to
116 * an integral value and adds that value to the exponent of the operand 'a'.
117 * The operation performed according to the IEC/IEEE Standard for Binary
118 * Floating-Point Arithmetic.
119 */
120
121 floatx80 floatx80_scale(floatx80 a, floatx80 b, float_status *status)
122 {
123 bool aSign, bSign;
124 int32_t aExp, bExp, shiftCount;
125 uint64_t aSig, bSig;
126
127 aSig = extractFloatx80Frac(a);
128 aExp = extractFloatx80Exp(a);
129 aSign = extractFloatx80Sign(a);
130 bSig = extractFloatx80Frac(b);
131 bExp = extractFloatx80Exp(b);
132 bSign = extractFloatx80Sign(b);
133
134 if (bExp == 0x7FFF) {
135 if ((uint64_t) (bSig << 1) ||
136 ((aExp == 0x7FFF) && (uint64_t) (aSig << 1))) {
137 return propagateFloatx80NaN(a, b, status);
138 }
139 float_raise(float_flag_invalid , status);
140 return floatx80_default_nan(status);
141 }
142 if (aExp == 0x7FFF) {
143 if ((uint64_t) (aSig << 1)) {
144 return propagateFloatx80NaN(a, b, status);
145 }
146 return floatx80_default_inf(aSign, status);
147 }
148 if (aExp == 0) {
149 if (aSig == 0) {
150 return packFloatx80(aSign, 0, 0);
151 }
152 if (bExp < 0x3FFF) {
153 return a;
154 }
155 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
156 }
157
158 if (bExp < 0x3FFF) {
159 return a;
160 }
161
162 if (0x400F < bExp) {
163 aExp = bSign ? -0x6001 : 0xE000;
164 return roundAndPackFloatx80(get_floatx80_rounding_precision(status),
165 aSign, aExp, aSig, 0, status);
166 }
167
168 shiftCount = 0x403E - bExp;
169 bSig >>= shiftCount;
170 aExp = bSign ? (aExp - bSig) : (aExp + bSig);
171
172 return roundAndPackFloatx80(get_floatx80_rounding_precision(status),
173 aSign, aExp, aSig, 0, status);
174 }
175
176 floatx80 floatx80_move(floatx80 a, float_status *status)
177 {
178 bool aSign;
179 int32_t aExp;
180 uint64_t aSig;
181
182 aSig = extractFloatx80Frac(a);
183 aExp = extractFloatx80Exp(a);
184 aSign = extractFloatx80Sign(a);
185
186 if (aExp == 0x7FFF) {
187 if ((uint64_t)(aSig << 1)) {
188 return propagateFloatx80NaNOneArg(a, status);
189 }
190 return a;
191 }
192 if (aExp == 0) {
193 if (aSig == 0) {
194 return a;
195 }
196 normalizeRoundAndPackFloatx80(get_floatx80_rounding_precision(status),
197 aSign, aExp, aSig, 0, status);
198 }
199 return roundAndPackFloatx80(get_floatx80_rounding_precision(status),
200 aSign,
201 aExp, aSig, 0, status);
202 }
203
204 /*
205 * Algorithms for transcendental functions supported by MC68881 and MC68882
206 * mathematical coprocessors. The functions are derived from FPSP library.
207 */
208
209 #define one_exp 0x3FFF
210 #define one_sig UINT64_C(0x8000000000000000)
211
212 /*
213 * Function for compactifying extended double-precision floating point values.
214 */
215
216 static int32_t floatx80_make_compact(int32_t aExp, uint64_t aSig)
217 {
218 return (aExp << 16) | (aSig >> 48);
219 }
220
221 /*
222 * Log base e of x plus 1
223 */
224
225 floatx80 floatx80_lognp1(floatx80 a, float_status *status)
226 {
227 bool aSign;
228 int32_t aExp;
229 uint64_t aSig, fSig;
230
231 FloatRoundMode user_rnd_mode;
232 FloatX80RoundPrec user_rnd_prec;
233
234 int32_t compact, j, k;
235 floatx80 fp0, fp1, fp2, fp3, f, logof2, klog2, saveu;
236
237 aSig = extractFloatx80Frac(a);
238 aExp = extractFloatx80Exp(a);
239 aSign = extractFloatx80Sign(a);
240
241 if (aExp == 0x7FFF) {
242 if ((uint64_t) (aSig << 1)) {
243 propagateFloatx80NaNOneArg(a, status);
244 }
245 if (aSign) {
246 float_raise(float_flag_invalid, status);
247 return floatx80_default_nan(status);
248 }
249 return floatx80_default_inf(0, status);
250 }
251
252 if (aExp == 0 && aSig == 0) {
253 return packFloatx80(aSign, 0, 0);
254 }
255
256 if (aSign && aExp >= one_exp) {
257 if (aExp == one_exp && aSig == one_sig) {
258 float_raise(float_flag_divbyzero, status);
259 return floatx80_default_inf(aSign, status);
260 }
261 float_raise(float_flag_invalid, status);
262 return floatx80_default_nan(status);
263 }
264
265 if (aExp < 0x3f99 || (aExp == 0x3f99 && aSig == one_sig)) {
266 /* <= min threshold */
267 float_raise(float_flag_inexact, status);
268 return floatx80_move(a, status);
269 }
270
271 user_rnd_mode = get_float_rounding_mode(status);
272 user_rnd_prec = get_floatx80_rounding_precision(status);
273 set_float_rounding_mode(float_round_nearest_even, status);
274 set_floatx80_rounding_precision(floatx80_precision_x, status);
275
276 compact = floatx80_make_compact(aExp, aSig);
277
278 fp0 = a; /* Z */
279 fp1 = a;
280
281 fp0 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000),
282 status), status); /* X = (1+Z) */
283
284 aExp = extractFloatx80Exp(fp0);
285 aSig = extractFloatx80Frac(fp0);
286
287 compact = floatx80_make_compact(aExp, aSig);
288
289 if (compact < 0x3FFE8000 || compact > 0x3FFFC000) {
290 /* |X| < 1/2 or |X| > 3/2 */
291 k = aExp - 0x3FFF;
292 fp1 = int32_to_floatx80(k, status);
293
294 fSig = (aSig & UINT64_C(0xFE00000000000000)) | UINT64_C(0x0100000000000000);
295 j = (fSig >> 56) & 0x7E; /* DISPLACEMENT FOR 1/F */
296
297 f = packFloatx80(0, 0x3FFF, fSig); /* F */
298 fp0 = packFloatx80(0, 0x3FFF, aSig); /* Y */
299
300 fp0 = floatx80_sub(fp0, f, status); /* Y-F */
301
302 lp1cont1:
303 /* LP1CONT1 */
304 fp0 = floatx80_mul(fp0, log_tbl[j], status); /* FP0 IS U = (Y-F)/F */
305 logof2 = packFloatx80(0, 0x3FFE, UINT64_C(0xB17217F7D1CF79AC));
306 klog2 = floatx80_mul(fp1, logof2, status); /* FP1 IS K*LOG2 */
307 fp2 = floatx80_mul(fp0, fp0, status); /* FP2 IS V=U*U */
308
309 fp3 = fp2;
310 fp1 = fp2;
311
312 fp1 = floatx80_mul(fp1, float64_to_floatx80(
313 make_float64(0x3FC2499AB5E4040B), status),
314 status); /* V*A6 */
315 fp2 = floatx80_mul(fp2, float64_to_floatx80(
316 make_float64(0xBFC555B5848CB7DB), status),
317 status); /* V*A5 */
318 fp1 = floatx80_add(fp1, float64_to_floatx80(
319 make_float64(0x3FC99999987D8730), status),
320 status); /* A4+V*A6 */
321 fp2 = floatx80_add(fp2, float64_to_floatx80(
322 make_float64(0xBFCFFFFFFF6F7E97), status),
323 status); /* A3+V*A5 */
324 fp1 = floatx80_mul(fp1, fp3, status); /* V*(A4+V*A6) */
325 fp2 = floatx80_mul(fp2, fp3, status); /* V*(A3+V*A5) */
326 fp1 = floatx80_add(fp1, float64_to_floatx80(
327 make_float64(0x3FD55555555555A4), status),
328 status); /* A2+V*(A4+V*A6) */
329 fp2 = floatx80_add(fp2, float64_to_floatx80(
330 make_float64(0xBFE0000000000008), status),
331 status); /* A1+V*(A3+V*A5) */
332 fp1 = floatx80_mul(fp1, fp3, status); /* V*(A2+V*(A4+V*A6)) */
333 fp2 = floatx80_mul(fp2, fp3, status); /* V*(A1+V*(A3+V*A5)) */
334 fp1 = floatx80_mul(fp1, fp0, status); /* U*V*(A2+V*(A4+V*A6)) */
335 fp0 = floatx80_add(fp0, fp2, status); /* U+V*(A1+V*(A3+V*A5)) */
336
337 fp1 = floatx80_add(fp1, log_tbl[j + 1],
338 status); /* LOG(F)+U*V*(A2+V*(A4+V*A6)) */
339 fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS LOG(F) + LOG(1+U) */
340
341 set_float_rounding_mode(user_rnd_mode, status);
342 set_floatx80_rounding_precision(user_rnd_prec, status);
343
344 a = floatx80_add(fp0, klog2, status);
345
346 float_raise(float_flag_inexact, status);
347
348 return a;
349 } else if (compact < 0x3FFEF07D || compact > 0x3FFF8841) {
350 /* |X| < 1/16 or |X| > -1/16 */
351 /* LP1CARE */
352 fSig = (aSig & UINT64_C(0xFE00000000000000)) | UINT64_C(0x0100000000000000);
353 f = packFloatx80(0, 0x3FFF, fSig); /* F */
354 j = (fSig >> 56) & 0x7E; /* DISPLACEMENT FOR 1/F */
355
356 if (compact >= 0x3FFF8000) { /* 1+Z >= 1 */
357 /* KISZERO */
358 fp0 = floatx80_sub(float32_to_floatx80(make_float32(0x3F800000),
359 status), f, status); /* 1-F */
360 fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS Y-F = (1-F)+Z */
361 fp1 = packFloatx80(0, 0, 0); /* K = 0 */
362 } else {
363 /* KISNEG */
364 fp0 = floatx80_sub(float32_to_floatx80(make_float32(0x40000000),
365 status), f, status); /* 2-F */
366 fp1 = floatx80_add(fp1, fp1, status); /* 2Z */
367 fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS Y-F = (2-F)+2Z */
368 fp1 = packFloatx80(1, one_exp, one_sig); /* K = -1 */
369 }
370 goto lp1cont1;
371 } else {
372 /* LP1ONE16 */
373 fp1 = floatx80_add(fp1, fp1, status); /* FP1 IS 2Z */
374 fp0 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000),
375 status), status); /* FP0 IS 1+X */
376
377 /* LP1CONT2 */
378 fp1 = floatx80_div(fp1, fp0, status); /* U */
379 saveu = fp1;
380 fp0 = floatx80_mul(fp1, fp1, status); /* FP0 IS V = U*U */
381 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS W = V*V */
382
383 fp3 = float64_to_floatx80(make_float64(0x3F175496ADD7DAD6),
384 status); /* B5 */
385 fp2 = float64_to_floatx80(make_float64(0x3F3C71C2FE80C7E0),
386 status); /* B4 */
387 fp3 = floatx80_mul(fp3, fp1, status); /* W*B5 */
388 fp2 = floatx80_mul(fp2, fp1, status); /* W*B4 */
389 fp3 = floatx80_add(fp3, float64_to_floatx80(
390 make_float64(0x3F624924928BCCFF), status),
391 status); /* B3+W*B5 */
392 fp2 = floatx80_add(fp2, float64_to_floatx80(
393 make_float64(0x3F899999999995EC), status),
394 status); /* B2+W*B4 */
395 fp1 = floatx80_mul(fp1, fp3, status); /* W*(B3+W*B5) */
396 fp2 = floatx80_mul(fp2, fp0, status); /* V*(B2+W*B4) */
397 fp1 = floatx80_add(fp1, float64_to_floatx80(
398 make_float64(0x3FB5555555555555), status),
399 status); /* B1+W*(B3+W*B5) */
400
401 fp0 = floatx80_mul(fp0, saveu, status); /* FP0 IS U*V */
402 fp1 = floatx80_add(fp1, fp2,
403 status); /* B1+W*(B3+W*B5) + V*(B2+W*B4) */
404 fp0 = floatx80_mul(fp0, fp1,
405 status); /* U*V*([B1+W*(B3+W*B5)] + [V*(B2+W*B4)]) */
406
407 set_float_rounding_mode(user_rnd_mode, status);
408 set_floatx80_rounding_precision(user_rnd_prec, status);
409
410 a = floatx80_add(fp0, saveu, status);
411
412 /*if (!floatx80_is_zero(a)) { */
413 float_raise(float_flag_inexact, status);
414 /*} */
415
416 return a;
417 }
418 }
419
420 /*
421 * Log base e
422 */
423
424 floatx80 floatx80_logn(floatx80 a, float_status *status)
425 {
426 bool aSign;
427 int32_t aExp;
428 uint64_t aSig, fSig;
429
430 FloatRoundMode user_rnd_mode;
431 FloatX80RoundPrec user_rnd_prec;
432
433 int32_t compact, j, k, adjk;
434 floatx80 fp0, fp1, fp2, fp3, f, logof2, klog2, saveu;
435
436 aSig = extractFloatx80Frac(a);
437 aExp = extractFloatx80Exp(a);
438 aSign = extractFloatx80Sign(a);
439
440 if (aExp == 0x7FFF) {
441 if ((uint64_t) (aSig << 1)) {
442 propagateFloatx80NaNOneArg(a, status);
443 }
444 if (aSign == 0) {
445 return floatx80_default_inf(0, status);
446 }
447 }
448
449 adjk = 0;
450
451 if (aExp == 0) {
452 if (aSig == 0) { /* zero */
453 float_raise(float_flag_divbyzero, status);
454 return floatx80_default_inf(1, status);
455 }
456 if ((aSig & one_sig) == 0) { /* denormal */
457 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
458 adjk = -100;
459 aExp += 100;
460 a = packFloatx80(aSign, aExp, aSig);
461 }
462 }
463
464 if (aSign) {
465 float_raise(float_flag_invalid, status);
466 return floatx80_default_nan(status);
467 }
468
469 user_rnd_mode = get_float_rounding_mode(status);
470 user_rnd_prec = get_floatx80_rounding_precision(status);
471 set_float_rounding_mode(float_round_nearest_even, status);
472 set_floatx80_rounding_precision(floatx80_precision_x, status);
473
474 compact = floatx80_make_compact(aExp, aSig);
475
476 if (compact < 0x3FFEF07D || compact > 0x3FFF8841) {
477 /* |X| < 15/16 or |X| > 17/16 */
478 k = aExp - 0x3FFF;
479 k += adjk;
480 fp1 = int32_to_floatx80(k, status);
481
482 fSig = (aSig & UINT64_C(0xFE00000000000000)) | UINT64_C(0x0100000000000000);
483 j = (fSig >> 56) & 0x7E; /* DISPLACEMENT FOR 1/F */
484
485 f = packFloatx80(0, 0x3FFF, fSig); /* F */
486 fp0 = packFloatx80(0, 0x3FFF, aSig); /* Y */
487
488 fp0 = floatx80_sub(fp0, f, status); /* Y-F */
489
490 /* LP1CONT1 */
491 fp0 = floatx80_mul(fp0, log_tbl[j], status); /* FP0 IS U = (Y-F)/F */
492 logof2 = packFloatx80(0, 0x3FFE, UINT64_C(0xB17217F7D1CF79AC));
493 klog2 = floatx80_mul(fp1, logof2, status); /* FP1 IS K*LOG2 */
494 fp2 = floatx80_mul(fp0, fp0, status); /* FP2 IS V=U*U */
495
496 fp3 = fp2;
497 fp1 = fp2;
498
499 fp1 = floatx80_mul(fp1, float64_to_floatx80(
500 make_float64(0x3FC2499AB5E4040B), status),
501 status); /* V*A6 */
502 fp2 = floatx80_mul(fp2, float64_to_floatx80(
503 make_float64(0xBFC555B5848CB7DB), status),
504 status); /* V*A5 */
505 fp1 = floatx80_add(fp1, float64_to_floatx80(
506 make_float64(0x3FC99999987D8730), status),
507 status); /* A4+V*A6 */
508 fp2 = floatx80_add(fp2, float64_to_floatx80(
509 make_float64(0xBFCFFFFFFF6F7E97), status),
510 status); /* A3+V*A5 */
511 fp1 = floatx80_mul(fp1, fp3, status); /* V*(A4+V*A6) */
512 fp2 = floatx80_mul(fp2, fp3, status); /* V*(A3+V*A5) */
513 fp1 = floatx80_add(fp1, float64_to_floatx80(
514 make_float64(0x3FD55555555555A4), status),
515 status); /* A2+V*(A4+V*A6) */
516 fp2 = floatx80_add(fp2, float64_to_floatx80(
517 make_float64(0xBFE0000000000008), status),
518 status); /* A1+V*(A3+V*A5) */
519 fp1 = floatx80_mul(fp1, fp3, status); /* V*(A2+V*(A4+V*A6)) */
520 fp2 = floatx80_mul(fp2, fp3, status); /* V*(A1+V*(A3+V*A5)) */
521 fp1 = floatx80_mul(fp1, fp0, status); /* U*V*(A2+V*(A4+V*A6)) */
522 fp0 = floatx80_add(fp0, fp2, status); /* U+V*(A1+V*(A3+V*A5)) */
523
524 fp1 = floatx80_add(fp1, log_tbl[j + 1],
525 status); /* LOG(F)+U*V*(A2+V*(A4+V*A6)) */
526 fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS LOG(F) + LOG(1+U) */
527
528 set_float_rounding_mode(user_rnd_mode, status);
529 set_floatx80_rounding_precision(user_rnd_prec, status);
530
531 a = floatx80_add(fp0, klog2, status);
532
533 float_raise(float_flag_inexact, status);
534
535 return a;
536 } else { /* |X-1| >= 1/16 */
537 fp0 = a;
538 fp1 = a;
539 fp1 = floatx80_sub(fp1, float32_to_floatx80(make_float32(0x3F800000),
540 status), status); /* FP1 IS X-1 */
541 fp0 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000),
542 status), status); /* FP0 IS X+1 */
543 fp1 = floatx80_add(fp1, fp1, status); /* FP1 IS 2(X-1) */
544
545 /* LP1CONT2 */
546 fp1 = floatx80_div(fp1, fp0, status); /* U */
547 saveu = fp1;
548 fp0 = floatx80_mul(fp1, fp1, status); /* FP0 IS V = U*U */
549 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS W = V*V */
550
551 fp3 = float64_to_floatx80(make_float64(0x3F175496ADD7DAD6),
552 status); /* B5 */
553 fp2 = float64_to_floatx80(make_float64(0x3F3C71C2FE80C7E0),
554 status); /* B4 */
555 fp3 = floatx80_mul(fp3, fp1, status); /* W*B5 */
556 fp2 = floatx80_mul(fp2, fp1, status); /* W*B4 */
557 fp3 = floatx80_add(fp3, float64_to_floatx80(
558 make_float64(0x3F624924928BCCFF), status),
559 status); /* B3+W*B5 */
560 fp2 = floatx80_add(fp2, float64_to_floatx80(
561 make_float64(0x3F899999999995EC), status),
562 status); /* B2+W*B4 */
563 fp1 = floatx80_mul(fp1, fp3, status); /* W*(B3+W*B5) */
564 fp2 = floatx80_mul(fp2, fp0, status); /* V*(B2+W*B4) */
565 fp1 = floatx80_add(fp1, float64_to_floatx80(
566 make_float64(0x3FB5555555555555), status),
567 status); /* B1+W*(B3+W*B5) */
568
569 fp0 = floatx80_mul(fp0, saveu, status); /* FP0 IS U*V */
570 fp1 = floatx80_add(fp1, fp2, status); /* B1+W*(B3+W*B5) + V*(B2+W*B4) */
571 fp0 = floatx80_mul(fp0, fp1,
572 status); /* U*V*([B1+W*(B3+W*B5)] + [V*(B2+W*B4)]) */
573
574 set_float_rounding_mode(user_rnd_mode, status);
575 set_floatx80_rounding_precision(user_rnd_prec, status);
576
577 a = floatx80_add(fp0, saveu, status);
578
579 /*if (!floatx80_is_zero(a)) { */
580 float_raise(float_flag_inexact, status);
581 /*} */
582
583 return a;
584 }
585 }
586
587 /*
588 * Log base 10
589 */
590
591 floatx80 floatx80_log10(floatx80 a, float_status *status)
592 {
593 bool aSign;
594 int32_t aExp;
595 uint64_t aSig;
596
597 FloatRoundMode user_rnd_mode;
598 FloatX80RoundPrec user_rnd_prec;
599
600 floatx80 fp0, fp1;
601
602 aSig = extractFloatx80Frac(a);
603 aExp = extractFloatx80Exp(a);
604 aSign = extractFloatx80Sign(a);
605
606 if (aExp == 0x7FFF) {
607 if ((uint64_t) (aSig << 1)) {
608 propagateFloatx80NaNOneArg(a, status);
609 }
610 if (aSign == 0) {
611 return floatx80_default_inf(0, status);
612 }
613 }
614
615 if (aExp == 0 && aSig == 0) {
616 float_raise(float_flag_divbyzero, status);
617 return floatx80_default_inf(1, status);
618 }
619
620 if (aSign) {
621 float_raise(float_flag_invalid, status);
622 return floatx80_default_nan(status);
623 }
624
625 user_rnd_mode = get_float_rounding_mode(status);
626 user_rnd_prec = get_floatx80_rounding_precision(status);
627 set_float_rounding_mode(float_round_nearest_even, status);
628 set_floatx80_rounding_precision(floatx80_precision_x, status);
629
630 fp0 = floatx80_logn(a, status);
631 fp1 = packFloatx80(0, 0x3FFD, UINT64_C(0xDE5BD8A937287195)); /* INV_L10 */
632
633 set_float_rounding_mode(user_rnd_mode, status);
634 set_floatx80_rounding_precision(user_rnd_prec, status);
635
636 a = floatx80_mul(fp0, fp1, status); /* LOGN(X)*INV_L10 */
637
638 float_raise(float_flag_inexact, status);
639
640 return a;
641 }
642
643 /*
644 * Log base 2
645 */
646
647 floatx80 floatx80_log2(floatx80 a, float_status *status)
648 {
649 bool aSign;
650 int32_t aExp;
651 uint64_t aSig;
652
653 FloatRoundMode user_rnd_mode;
654 FloatX80RoundPrec user_rnd_prec;
655
656 floatx80 fp0, fp1;
657
658 aSig = extractFloatx80Frac(a);
659 aExp = extractFloatx80Exp(a);
660 aSign = extractFloatx80Sign(a);
661
662 if (aExp == 0x7FFF) {
663 if ((uint64_t) (aSig << 1)) {
664 propagateFloatx80NaNOneArg(a, status);
665 }
666 if (aSign == 0) {
667 return floatx80_default_inf(0, status);
668 }
669 }
670
671 if (aExp == 0) {
672 if (aSig == 0) {
673 float_raise(float_flag_divbyzero, status);
674 return floatx80_default_inf(1, status);
675 }
676 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
677 }
678
679 if (aSign) {
680 float_raise(float_flag_invalid, status);
681 return floatx80_default_nan(status);
682 }
683
684 user_rnd_mode = get_float_rounding_mode(status);
685 user_rnd_prec = get_floatx80_rounding_precision(status);
686 set_float_rounding_mode(float_round_nearest_even, status);
687 set_floatx80_rounding_precision(floatx80_precision_x, status);
688
689 if (aSig == one_sig) { /* X is 2^k */
690 set_float_rounding_mode(user_rnd_mode, status);
691 set_floatx80_rounding_precision(user_rnd_prec, status);
692
693 a = int32_to_floatx80(aExp - 0x3FFF, status);
694 } else {
695 fp0 = floatx80_logn(a, status);
696 fp1 = packFloatx80(0, 0x3FFF, UINT64_C(0xB8AA3B295C17F0BC)); /* INV_L2 */
697
698 set_float_rounding_mode(user_rnd_mode, status);
699 set_floatx80_rounding_precision(user_rnd_prec, status);
700
701 a = floatx80_mul(fp0, fp1, status); /* LOGN(X)*INV_L2 */
702 }
703
704 float_raise(float_flag_inexact, status);
705
706 return a;
707 }
708
709 /*
710 * e to x
711 */
712
713 floatx80 floatx80_etox(floatx80 a, float_status *status)
714 {
715 bool aSign;
716 int32_t aExp;
717 uint64_t aSig;
718
719 FloatRoundMode user_rnd_mode;
720 FloatX80RoundPrec user_rnd_prec;
721
722 int32_t compact, n, j, k, m, m1;
723 floatx80 fp0, fp1, fp2, fp3, l2, scale, adjscale;
724 bool adjflag;
725
726 aSig = extractFloatx80Frac(a);
727 aExp = extractFloatx80Exp(a);
728 aSign = extractFloatx80Sign(a);
729
730 if (aExp == 0x7FFF) {
731 if ((uint64_t) (aSig << 1)) {
732 return propagateFloatx80NaNOneArg(a, status);
733 }
734 if (aSign) {
735 return packFloatx80(0, 0, 0);
736 }
737 return floatx80_default_inf(0, status);
738 }
739
740 if (aExp == 0 && aSig == 0) {
741 return packFloatx80(0, one_exp, one_sig);
742 }
743
744 user_rnd_mode = get_float_rounding_mode(status);
745 user_rnd_prec = get_floatx80_rounding_precision(status);
746 set_float_rounding_mode(float_round_nearest_even, status);
747 set_floatx80_rounding_precision(floatx80_precision_x, status);
748
749 adjflag = 0;
750
751 if (aExp >= 0x3FBE) { /* |X| >= 2^(-65) */
752 compact = floatx80_make_compact(aExp, aSig);
753
754 if (compact < 0x400CB167) { /* |X| < 16380 log2 */
755 fp0 = a;
756 fp1 = a;
757 fp0 = floatx80_mul(fp0, float32_to_floatx80(
758 make_float32(0x42B8AA3B), status),
759 status); /* 64/log2 * X */
760 adjflag = 0;
761 n = floatx80_to_int32(fp0, status); /* int(64/log2*X) */
762 fp0 = int32_to_floatx80(n, status);
763
764 j = n & 0x3F; /* J = N mod 64 */
765 m = n / 64; /* NOTE: this is really arithmetic right shift by 6 */
766 if (n < 0 && j) {
767 /*
768 * arithmetic right shift is division and
769 * round towards minus infinity
770 */
771 m--;
772 }
773 m += 0x3FFF; /* biased exponent of 2^(M) */
774
775 expcont1:
776 fp2 = fp0; /* N */
777 fp0 = floatx80_mul(fp0, float32_to_floatx80(
778 make_float32(0xBC317218), status),
779 status); /* N * L1, L1 = lead(-log2/64) */
780 l2 = packFloatx80(0, 0x3FDC, UINT64_C(0x82E308654361C4C6));
781 fp2 = floatx80_mul(fp2, l2, status); /* N * L2, L1+L2 = -log2/64 */
782 fp0 = floatx80_add(fp0, fp1, status); /* X + N*L1 */
783 fp0 = floatx80_add(fp0, fp2, status); /* R */
784
785 fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */
786 fp2 = float32_to_floatx80(make_float32(0x3AB60B70),
787 status); /* A5 */
788 fp2 = floatx80_mul(fp2, fp1, status); /* fp2 is S*A5 */
789 fp3 = floatx80_mul(float32_to_floatx80(make_float32(0x3C088895),
790 status), fp1,
791 status); /* fp3 is S*A4 */
792 fp2 = floatx80_add(fp2, float64_to_floatx80(make_float64(
793 0x3FA5555555554431), status),
794 status); /* fp2 is A3+S*A5 */
795 fp3 = floatx80_add(fp3, float64_to_floatx80(make_float64(
796 0x3FC5555555554018), status),
797 status); /* fp3 is A2+S*A4 */
798 fp2 = floatx80_mul(fp2, fp1, status); /* fp2 is S*(A3+S*A5) */
799 fp3 = floatx80_mul(fp3, fp1, status); /* fp3 is S*(A2+S*A4) */
800 fp2 = floatx80_add(fp2, float32_to_floatx80(
801 make_float32(0x3F000000), status),
802 status); /* fp2 is A1+S*(A3+S*A5) */
803 fp3 = floatx80_mul(fp3, fp0, status); /* fp3 IS R*S*(A2+S*A4) */
804 fp2 = floatx80_mul(fp2, fp1,
805 status); /* fp2 IS S*(A1+S*(A3+S*A5)) */
806 fp0 = floatx80_add(fp0, fp3, status); /* fp0 IS R+R*S*(A2+S*A4) */
807 fp0 = floatx80_add(fp0, fp2, status); /* fp0 IS EXP(R) - 1 */
808
809 fp1 = exp_tbl[j];
810 fp0 = floatx80_mul(fp0, fp1, status); /* 2^(J/64)*(Exp(R)-1) */
811 fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j], status),
812 status); /* accurate 2^(J/64) */
813 fp0 = floatx80_add(fp0, fp1,
814 status); /* 2^(J/64) + 2^(J/64)*(Exp(R)-1) */
815
816 scale = packFloatx80(0, m, one_sig);
817 if (adjflag) {
818 adjscale = packFloatx80(0, m1, one_sig);
819 fp0 = floatx80_mul(fp0, adjscale, status);
820 }
821
822 set_float_rounding_mode(user_rnd_mode, status);
823 set_floatx80_rounding_precision(user_rnd_prec, status);
824
825 a = floatx80_mul(fp0, scale, status);
826
827 float_raise(float_flag_inexact, status);
828
829 return a;
830 } else { /* |X| >= 16380 log2 */
831 if (compact > 0x400CB27C) { /* |X| >= 16480 log2 */
832 set_float_rounding_mode(user_rnd_mode, status);
833 set_floatx80_rounding_precision(user_rnd_prec, status);
834 if (aSign) {
835 a = roundAndPackFloatx80(get_floatx80_rounding_precision(status),
836 0, -0x1000, aSig, 0, status);
837 } else {
838 a = roundAndPackFloatx80(get_floatx80_rounding_precision(status),
839 0, 0x8000, aSig, 0, status);
840 }
841 float_raise(float_flag_inexact, status);
842
843 return a;
844 } else {
845 fp0 = a;
846 fp1 = a;
847 fp0 = floatx80_mul(fp0, float32_to_floatx80(
848 make_float32(0x42B8AA3B), status),
849 status); /* 64/log2 * X */
850 adjflag = 1;
851 n = floatx80_to_int32(fp0, status); /* int(64/log2*X) */
852 fp0 = int32_to_floatx80(n, status);
853
854 j = n & 0x3F; /* J = N mod 64 */
855 /* NOTE: this is really arithmetic right shift by 6 */
856 k = n / 64;
857 if (n < 0 && j) {
858 /* arithmetic right shift is division and
859 * round towards minus infinity
860 */
861 k--;
862 }
863 /* NOTE: this is really arithmetic right shift by 1 */
864 m1 = k / 2;
865 if (k < 0 && (k & 1)) {
866 /* arithmetic right shift is division and
867 * round towards minus infinity
868 */
869 m1--;
870 }
871 m = k - m1;
872 m1 += 0x3FFF; /* biased exponent of 2^(M1) */
873 m += 0x3FFF; /* biased exponent of 2^(M) */
874
875 goto expcont1;
876 }
877 }
878 } else { /* |X| < 2^(-65) */
879 set_float_rounding_mode(user_rnd_mode, status);
880 set_floatx80_rounding_precision(user_rnd_prec, status);
881
882 a = floatx80_add(a, float32_to_floatx80(make_float32(0x3F800000),
883 status), status); /* 1 + X */
884
885 float_raise(float_flag_inexact, status);
886
887 return a;
888 }
889 }
890
891 /*
892 * 2 to x
893 */
894
895 floatx80 floatx80_twotox(floatx80 a, float_status *status)
896 {
897 bool aSign;
898 int32_t aExp;
899 uint64_t aSig;
900
901 FloatRoundMode user_rnd_mode;
902 FloatX80RoundPrec user_rnd_prec;
903
904 int32_t compact, n, j, l, m, m1;
905 floatx80 fp0, fp1, fp2, fp3, adjfact, fact1, fact2;
906
907 aSig = extractFloatx80Frac(a);
908 aExp = extractFloatx80Exp(a);
909 aSign = extractFloatx80Sign(a);
910
911 if (aExp == 0x7FFF) {
912 if ((uint64_t) (aSig << 1)) {
913 return propagateFloatx80NaNOneArg(a, status);
914 }
915 if (aSign) {
916 return packFloatx80(0, 0, 0);
917 }
918 return floatx80_default_inf(0, status);
919 }
920
921 if (aExp == 0 && aSig == 0) {
922 return packFloatx80(0, one_exp, one_sig);
923 }
924
925 user_rnd_mode = get_float_rounding_mode(status);
926 user_rnd_prec = get_floatx80_rounding_precision(status);
927 set_float_rounding_mode(float_round_nearest_even, status);
928 set_floatx80_rounding_precision(floatx80_precision_x, status);
929
930 fp0 = a;
931
932 compact = floatx80_make_compact(aExp, aSig);
933
934 if (compact < 0x3FB98000 || compact > 0x400D80C0) {
935 /* |X| > 16480 or |X| < 2^(-70) */
936 if (compact > 0x3FFF8000) { /* |X| > 16480 */
937 set_float_rounding_mode(user_rnd_mode, status);
938 set_floatx80_rounding_precision(user_rnd_prec, status);
939
940 if (aSign) {
941 return roundAndPackFloatx80(get_floatx80_rounding_precision(status),
942 0, -0x1000, aSig, 0, status);
943 } else {
944 return roundAndPackFloatx80(get_floatx80_rounding_precision(status),
945 0, 0x8000, aSig, 0, status);
946 }
947 } else { /* |X| < 2^(-70) */
948 set_float_rounding_mode(user_rnd_mode, status);
949 set_floatx80_rounding_precision(user_rnd_prec, status);
950
951 a = floatx80_add(fp0, float32_to_floatx80(
952 make_float32(0x3F800000), status),
953 status); /* 1 + X */
954
955 float_raise(float_flag_inexact, status);
956
957 return a;
958 }
959 } else { /* 2^(-70) <= |X| <= 16480 */
960 fp1 = fp0; /* X */
961 fp1 = floatx80_mul(fp1, float32_to_floatx80(
962 make_float32(0x42800000), status),
963 status); /* X * 64 */
964 n = floatx80_to_int32(fp1, status);
965 fp1 = int32_to_floatx80(n, status);
966 j = n & 0x3F;
967 l = n / 64; /* NOTE: this is really arithmetic right shift by 6 */
968 if (n < 0 && j) {
969 /*
970 * arithmetic right shift is division and
971 * round towards minus infinity
972 */
973 l--;
974 }
975 m = l / 2; /* NOTE: this is really arithmetic right shift by 1 */
976 if (l < 0 && (l & 1)) {
977 /*
978 * arithmetic right shift is division and
979 * round towards minus infinity
980 */
981 m--;
982 }
983 m1 = l - m;
984 m1 += 0x3FFF; /* ADJFACT IS 2^(M') */
985
986 adjfact = packFloatx80(0, m1, one_sig);
987 fact1 = exp2_tbl[j];
988 fact1.high += m;
989 fact2.high = exp2_tbl2[j] >> 16;
990 fact2.high += m;
991 fact2.low = (uint64_t)(exp2_tbl2[j] & 0xFFFF);
992 fact2.low <<= 48;
993
994 fp1 = floatx80_mul(fp1, float32_to_floatx80(
995 make_float32(0x3C800000), status),
996 status); /* (1/64)*N */
997 fp0 = floatx80_sub(fp0, fp1, status); /* X - (1/64)*INT(64 X) */
998 fp2 = packFloatx80(0, 0x3FFE, UINT64_C(0xB17217F7D1CF79AC)); /* LOG2 */
999 fp0 = floatx80_mul(fp0, fp2, status); /* R */
1000
1001 /* EXPR */
1002 fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */
1003 fp2 = float64_to_floatx80(make_float64(0x3F56C16D6F7BD0B2),
1004 status); /* A5 */
1005 fp3 = float64_to_floatx80(make_float64(0x3F811112302C712C),
1006 status); /* A4 */
1007 fp2 = floatx80_mul(fp2, fp1, status); /* S*A5 */
1008 fp3 = floatx80_mul(fp3, fp1, status); /* S*A4 */
1009 fp2 = floatx80_add(fp2, float64_to_floatx80(
1010 make_float64(0x3FA5555555554CC1), status),
1011 status); /* A3+S*A5 */
1012 fp3 = floatx80_add(fp3, float64_to_floatx80(
1013 make_float64(0x3FC5555555554A54), status),
1014 status); /* A2+S*A4 */
1015 fp2 = floatx80_mul(fp2, fp1, status); /* S*(A3+S*A5) */
1016 fp3 = floatx80_mul(fp3, fp1, status); /* S*(A2+S*A4) */
1017 fp2 = floatx80_add(fp2, float64_to_floatx80(
1018 make_float64(0x3FE0000000000000), status),
1019 status); /* A1+S*(A3+S*A5) */
1020 fp3 = floatx80_mul(fp3, fp0, status); /* R*S*(A2+S*A4) */
1021
1022 fp2 = floatx80_mul(fp2, fp1, status); /* S*(A1+S*(A3+S*A5)) */
1023 fp0 = floatx80_add(fp0, fp3, status); /* R+R*S*(A2+S*A4) */
1024 fp0 = floatx80_add(fp0, fp2, status); /* EXP(R) - 1 */
1025
1026 fp0 = floatx80_mul(fp0, fact1, status);
1027 fp0 = floatx80_add(fp0, fact2, status);
1028 fp0 = floatx80_add(fp0, fact1, status);
1029
1030 set_float_rounding_mode(user_rnd_mode, status);
1031 set_floatx80_rounding_precision(user_rnd_prec, status);
1032
1033 a = floatx80_mul(fp0, adjfact, status);
1034
1035 float_raise(float_flag_inexact, status);
1036
1037 return a;
1038 }
1039 }
1040
1041 /*
1042 * 10 to x
1043 */
1044
1045 floatx80 floatx80_tentox(floatx80 a, float_status *status)
1046 {
1047 bool aSign;
1048 int32_t aExp;
1049 uint64_t aSig;
1050
1051 FloatRoundMode user_rnd_mode;
1052 FloatX80RoundPrec user_rnd_prec;
1053
1054 int32_t compact, n, j, l, m, m1;
1055 floatx80 fp0, fp1, fp2, fp3, adjfact, fact1, fact2;
1056
1057 aSig = extractFloatx80Frac(a);
1058 aExp = extractFloatx80Exp(a);
1059 aSign = extractFloatx80Sign(a);
1060
1061 if (aExp == 0x7FFF) {
1062 if ((uint64_t) (aSig << 1)) {
1063 return propagateFloatx80NaNOneArg(a, status);
1064 }
1065 if (aSign) {
1066 return packFloatx80(0, 0, 0);
1067 }
1068 return floatx80_default_inf(0, status);
1069 }
1070
1071 if (aExp == 0 && aSig == 0) {
1072 return packFloatx80(0, one_exp, one_sig);
1073 }
1074
1075 user_rnd_mode = get_float_rounding_mode(status);
1076 user_rnd_prec = get_floatx80_rounding_precision(status);
1077 set_float_rounding_mode(float_round_nearest_even, status);
1078 set_floatx80_rounding_precision(floatx80_precision_x, status);
1079
1080 fp0 = a;
1081
1082 compact = floatx80_make_compact(aExp, aSig);
1083
1084 if (compact < 0x3FB98000 || compact > 0x400B9B07) {
1085 /* |X| > 16480 LOG2/LOG10 or |X| < 2^(-70) */
1086 if (compact > 0x3FFF8000) { /* |X| > 16480 */
1087 set_float_rounding_mode(user_rnd_mode, status);
1088 set_floatx80_rounding_precision(user_rnd_prec, status);
1089
1090 if (aSign) {
1091 return roundAndPackFloatx80(get_floatx80_rounding_precision(status),
1092 0, -0x1000, aSig, 0, status);
1093 } else {
1094 return roundAndPackFloatx80(get_floatx80_rounding_precision(status),
1095 0, 0x8000, aSig, 0, status);
1096 }
1097 } else { /* |X| < 2^(-70) */
1098 set_float_rounding_mode(user_rnd_mode, status);
1099 set_floatx80_rounding_precision(user_rnd_prec, status);
1100
1101 a = floatx80_add(fp0, float32_to_floatx80(
1102 make_float32(0x3F800000), status),
1103 status); /* 1 + X */
1104
1105 float_raise(float_flag_inexact, status);
1106
1107 return a;
1108 }
1109 } else { /* 2^(-70) <= |X| <= 16480 LOG 2 / LOG 10 */
1110 fp1 = fp0; /* X */
1111 fp1 = floatx80_mul(fp1, float64_to_floatx80(
1112 make_float64(0x406A934F0979A371),
1113 status), status); /* X*64*LOG10/LOG2 */
1114 n = floatx80_to_int32(fp1, status); /* N=INT(X*64*LOG10/LOG2) */
1115 fp1 = int32_to_floatx80(n, status);
1116
1117 j = n & 0x3F;
1118 l = n / 64; /* NOTE: this is really arithmetic right shift by 6 */
1119 if (n < 0 && j) {
1120 /*
1121 * arithmetic right shift is division and
1122 * round towards minus infinity
1123 */
1124 l--;
1125 }
1126 m = l / 2; /* NOTE: this is really arithmetic right shift by 1 */
1127 if (l < 0 && (l & 1)) {
1128 /*
1129 * arithmetic right shift is division and
1130 * round towards minus infinity
1131 */
1132 m--;
1133 }
1134 m1 = l - m;
1135 m1 += 0x3FFF; /* ADJFACT IS 2^(M') */
1136
1137 adjfact = packFloatx80(0, m1, one_sig);
1138 fact1 = exp2_tbl[j];
1139 fact1.high += m;
1140 fact2.high = exp2_tbl2[j] >> 16;
1141 fact2.high += m;
1142 fact2.low = (uint64_t)(exp2_tbl2[j] & 0xFFFF);
1143 fact2.low <<= 48;
1144
1145 fp2 = fp1; /* N */
1146 fp1 = floatx80_mul(fp1, float64_to_floatx80(
1147 make_float64(0x3F734413509F8000), status),
1148 status); /* N*(LOG2/64LOG10)_LEAD */
1149 fp3 = packFloatx80(1, 0x3FCD, UINT64_C(0xC0219DC1DA994FD2));
1150 fp2 = floatx80_mul(fp2, fp3, status); /* N*(LOG2/64LOG10)_TRAIL */
1151 fp0 = floatx80_sub(fp0, fp1, status); /* X - N L_LEAD */
1152 fp0 = floatx80_sub(fp0, fp2, status); /* X - N L_TRAIL */
1153 fp2 = packFloatx80(0, 0x4000, UINT64_C(0x935D8DDDAAA8AC17)); /* LOG10 */
1154 fp0 = floatx80_mul(fp0, fp2, status); /* R */
1155
1156 /* EXPR */
1157 fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */
1158 fp2 = float64_to_floatx80(make_float64(0x3F56C16D6F7BD0B2),
1159 status); /* A5 */
1160 fp3 = float64_to_floatx80(make_float64(0x3F811112302C712C),
1161 status); /* A4 */
1162 fp2 = floatx80_mul(fp2, fp1, status); /* S*A5 */
1163 fp3 = floatx80_mul(fp3, fp1, status); /* S*A4 */
1164 fp2 = floatx80_add(fp2, float64_to_floatx80(
1165 make_float64(0x3FA5555555554CC1), status),
1166 status); /* A3+S*A5 */
1167 fp3 = floatx80_add(fp3, float64_to_floatx80(
1168 make_float64(0x3FC5555555554A54), status),
1169 status); /* A2+S*A4 */
1170 fp2 = floatx80_mul(fp2, fp1, status); /* S*(A3+S*A5) */
1171 fp3 = floatx80_mul(fp3, fp1, status); /* S*(A2+S*A4) */
1172 fp2 = floatx80_add(fp2, float64_to_floatx80(
1173 make_float64(0x3FE0000000000000), status),
1174 status); /* A1+S*(A3+S*A5) */
1175 fp3 = floatx80_mul(fp3, fp0, status); /* R*S*(A2+S*A4) */
1176
1177 fp2 = floatx80_mul(fp2, fp1, status); /* S*(A1+S*(A3+S*A5)) */
1178 fp0 = floatx80_add(fp0, fp3, status); /* R+R*S*(A2+S*A4) */
1179 fp0 = floatx80_add(fp0, fp2, status); /* EXP(R) - 1 */
1180
1181 fp0 = floatx80_mul(fp0, fact1, status);
1182 fp0 = floatx80_add(fp0, fact2, status);
1183 fp0 = floatx80_add(fp0, fact1, status);
1184
1185 set_float_rounding_mode(user_rnd_mode, status);
1186 set_floatx80_rounding_precision(user_rnd_prec, status);
1187
1188 a = floatx80_mul(fp0, adjfact, status);
1189
1190 float_raise(float_flag_inexact, status);
1191
1192 return a;
1193 }
1194 }
1195
1196 /*
1197 * Tangent
1198 */
1199
1200 floatx80 floatx80_tan(floatx80 a, float_status *status)
1201 {
1202 bool aSign, xSign;
1203 int32_t aExp, xExp;
1204 uint64_t aSig, xSig;
1205
1206 FloatRoundMode user_rnd_mode;
1207 FloatX80RoundPrec user_rnd_prec;
1208
1209 int32_t compact, l, n, j;
1210 floatx80 fp0, fp1, fp2, fp3, fp4, fp5, invtwopi, twopi1, twopi2;
1211 float32 twoto63;
1212 bool endflag;
1213
1214 aSig = extractFloatx80Frac(a);
1215 aExp = extractFloatx80Exp(a);
1216 aSign = extractFloatx80Sign(a);
1217
1218 if (aExp == 0x7FFF) {
1219 if ((uint64_t) (aSig << 1)) {
1220 return propagateFloatx80NaNOneArg(a, status);
1221 }
1222 float_raise(float_flag_invalid, status);
1223 return floatx80_default_nan(status);
1224 }
1225
1226 if (aExp == 0 && aSig == 0) {
1227 return packFloatx80(aSign, 0, 0);
1228 }
1229
1230 user_rnd_mode = get_float_rounding_mode(status);
1231 user_rnd_prec = get_floatx80_rounding_precision(status);
1232 set_float_rounding_mode(float_round_nearest_even, status);
1233 set_floatx80_rounding_precision(floatx80_precision_x, status);
1234
1235 compact = floatx80_make_compact(aExp, aSig);
1236
1237 fp0 = a;
1238
1239 if (compact < 0x3FD78000 || compact > 0x4004BC7E) {
1240 /* 2^(-40) > |X| > 15 PI */
1241 if (compact > 0x3FFF8000) { /* |X| >= 15 PI */
1242 /* REDUCEX */
1243 fp1 = packFloatx80(0, 0, 0);
1244 if (compact == 0x7FFEFFFF) {
1245 twopi1 = packFloatx80(aSign ^ 1, 0x7FFE,
1246 UINT64_C(0xC90FDAA200000000));
1247 twopi2 = packFloatx80(aSign ^ 1, 0x7FDC,
1248 UINT64_C(0x85A308D300000000));
1249 fp0 = floatx80_add(fp0, twopi1, status);
1250 fp1 = fp0;
1251 fp0 = floatx80_add(fp0, twopi2, status);
1252 fp1 = floatx80_sub(fp1, fp0, status);
1253 fp1 = floatx80_add(fp1, twopi2, status);
1254 }
1255 loop:
1256 xSign = extractFloatx80Sign(fp0);
1257 xExp = extractFloatx80Exp(fp0);
1258 xExp -= 0x3FFF;
1259 if (xExp <= 28) {
1260 l = 0;
1261 endflag = true;
1262 } else {
1263 l = xExp - 27;
1264 endflag = false;
1265 }
1266 invtwopi = packFloatx80(0, 0x3FFE - l,
1267 UINT64_C(0xA2F9836E4E44152A)); /* INVTWOPI */
1268 twopi1 = packFloatx80(0, 0x3FFF + l, UINT64_C(0xC90FDAA200000000));
1269 twopi2 = packFloatx80(0, 0x3FDD + l, UINT64_C(0x85A308D300000000));
1270
1271 /* SIGN(INARG)*2^63 IN SGL */
1272 twoto63 = packFloat32(xSign, 0xBE, 0);
1273
1274 fp2 = floatx80_mul(fp0, invtwopi, status);
1275 fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63, status),
1276 status); /* THE FRACT PART OF FP2 IS ROUNDED */
1277 fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63, status),
1278 status); /* FP2 is N */
1279 fp4 = floatx80_mul(twopi1, fp2, status); /* W = N*P1 */
1280 fp5 = floatx80_mul(twopi2, fp2, status); /* w = N*P2 */
1281 fp3 = floatx80_add(fp4, fp5, status); /* FP3 is P */
1282 fp4 = floatx80_sub(fp4, fp3, status); /* W-P */
1283 fp0 = floatx80_sub(fp0, fp3, status); /* FP0 is A := R - P */
1284 fp4 = floatx80_add(fp4, fp5, status); /* FP4 is p = (W-P)+w */
1285 fp3 = fp0; /* FP3 is A */
1286 fp1 = floatx80_sub(fp1, fp4, status); /* FP1 is a := r - p */
1287 fp0 = floatx80_add(fp0, fp1, status); /* FP0 is R := A+a */
1288
1289 if (endflag) {
1290 n = floatx80_to_int32(fp2, status);
1291 goto tancont;
1292 }
1293 fp3 = floatx80_sub(fp3, fp0, status); /* A-R */
1294 fp1 = floatx80_add(fp1, fp3, status); /* FP1 is r := (A-R)+a */
1295 goto loop;
1296 } else {
1297 set_float_rounding_mode(user_rnd_mode, status);
1298 set_floatx80_rounding_precision(user_rnd_prec, status);
1299
1300 a = floatx80_move(a, status);
1301
1302 float_raise(float_flag_inexact, status);
1303
1304 return a;
1305 }
1306 } else {
1307 fp1 = floatx80_mul(fp0, float64_to_floatx80(
1308 make_float64(0x3FE45F306DC9C883), status),
1309 status); /* X*2/PI */
1310
1311 n = floatx80_to_int32(fp1, status);
1312 j = 32 + n;
1313
1314 fp0 = floatx80_sub(fp0, pi_tbl[j], status); /* X-Y1 */
1315 fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j], status),
1316 status); /* FP0 IS R = (X-Y1)-Y2 */
1317
1318 tancont:
1319 if (n & 1) {
1320 /* NODD */
1321 fp1 = fp0; /* R */
1322 fp0 = floatx80_mul(fp0, fp0, status); /* S = R*R */
1323 fp3 = float64_to_floatx80(make_float64(0x3EA0B759F50F8688),
1324 status); /* Q4 */
1325 fp2 = float64_to_floatx80(make_float64(0xBEF2BAA5A8924F04),
1326 status); /* P3 */
1327 fp3 = floatx80_mul(fp3, fp0, status); /* SQ4 */
1328 fp2 = floatx80_mul(fp2, fp0, status); /* SP3 */
1329 fp3 = floatx80_add(fp3, float64_to_floatx80(
1330 make_float64(0xBF346F59B39BA65F), status),
1331 status); /* Q3+SQ4 */
1332 fp4 = packFloatx80(0, 0x3FF6, UINT64_C(0xE073D3FC199C4A00));
1333 fp2 = floatx80_add(fp2, fp4, status); /* P2+SP3 */
1334 fp3 = floatx80_mul(fp3, fp0, status); /* S(Q3+SQ4) */
1335 fp2 = floatx80_mul(fp2, fp0, status); /* S(P2+SP3) */
1336 fp4 = packFloatx80(0, 0x3FF9, UINT64_C(0xD23CD68415D95FA1));
1337 fp3 = floatx80_add(fp3, fp4, status); /* Q2+S(Q3+SQ4) */
1338 fp4 = packFloatx80(1, 0x3FFC, UINT64_C(0x8895A6C5FB423BCA));
1339 fp2 = floatx80_add(fp2, fp4, status); /* P1+S(P2+SP3) */
1340 fp3 = floatx80_mul(fp3, fp0, status); /* S(Q2+S(Q3+SQ4)) */
1341 fp2 = floatx80_mul(fp2, fp0, status); /* S(P1+S(P2+SP3)) */
1342 fp4 = packFloatx80(1, 0x3FFD, UINT64_C(0xEEF57E0DA84BC8CE));
1343 fp3 = floatx80_add(fp3, fp4, status); /* Q1+S(Q2+S(Q3+SQ4)) */
1344 fp2 = floatx80_mul(fp2, fp1, status); /* RS(P1+S(P2+SP3)) */
1345 fp0 = floatx80_mul(fp0, fp3, status); /* S(Q1+S(Q2+S(Q3+SQ4))) */
1346 fp1 = floatx80_add(fp1, fp2, status); /* R+RS(P1+S(P2+SP3)) */
1347 fp0 = floatx80_add(fp0, float32_to_floatx80(
1348 make_float32(0x3F800000), status),
1349 status); /* 1+S(Q1+S(Q2+S(Q3+SQ4))) */
1350
1351 xSign = extractFloatx80Sign(fp1);
1352 xExp = extractFloatx80Exp(fp1);
1353 xSig = extractFloatx80Frac(fp1);
1354 xSign ^= 1;
1355 fp1 = packFloatx80(xSign, xExp, xSig);
1356
1357 set_float_rounding_mode(user_rnd_mode, status);
1358 set_floatx80_rounding_precision(user_rnd_prec, status);
1359
1360 a = floatx80_div(fp0, fp1, status);
1361
1362 float_raise(float_flag_inexact, status);
1363
1364 return a;
1365 } else {
1366 fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */
1367 fp3 = float64_to_floatx80(make_float64(0x3EA0B759F50F8688),
1368 status); /* Q4 */
1369 fp2 = float64_to_floatx80(make_float64(0xBEF2BAA5A8924F04),
1370 status); /* P3 */
1371 fp3 = floatx80_mul(fp3, fp1, status); /* SQ4 */
1372 fp2 = floatx80_mul(fp2, fp1, status); /* SP3 */
1373 fp3 = floatx80_add(fp3, float64_to_floatx80(
1374 make_float64(0xBF346F59B39BA65F), status),
1375 status); /* Q3+SQ4 */
1376 fp4 = packFloatx80(0, 0x3FF6, UINT64_C(0xE073D3FC199C4A00));
1377 fp2 = floatx80_add(fp2, fp4, status); /* P2+SP3 */
1378 fp3 = floatx80_mul(fp3, fp1, status); /* S(Q3+SQ4) */
1379 fp2 = floatx80_mul(fp2, fp1, status); /* S(P2+SP3) */
1380 fp4 = packFloatx80(0, 0x3FF9, UINT64_C(0xD23CD68415D95FA1));
1381 fp3 = floatx80_add(fp3, fp4, status); /* Q2+S(Q3+SQ4) */
1382 fp4 = packFloatx80(1, 0x3FFC, UINT64_C(0x8895A6C5FB423BCA));
1383 fp2 = floatx80_add(fp2, fp4, status); /* P1+S(P2+SP3) */
1384 fp3 = floatx80_mul(fp3, fp1, status); /* S(Q2+S(Q3+SQ4)) */
1385 fp2 = floatx80_mul(fp2, fp1, status); /* S(P1+S(P2+SP3)) */
1386 fp4 = packFloatx80(1, 0x3FFD, UINT64_C(0xEEF57E0DA84BC8CE));
1387 fp3 = floatx80_add(fp3, fp4, status); /* Q1+S(Q2+S(Q3+SQ4)) */
1388 fp2 = floatx80_mul(fp2, fp0, status); /* RS(P1+S(P2+SP3)) */
1389 fp1 = floatx80_mul(fp1, fp3, status); /* S(Q1+S(Q2+S(Q3+SQ4))) */
1390 fp0 = floatx80_add(fp0, fp2, status); /* R+RS(P1+S(P2+SP3)) */
1391 fp1 = floatx80_add(fp1, float32_to_floatx80(
1392 make_float32(0x3F800000), status),
1393 status); /* 1+S(Q1+S(Q2+S(Q3+SQ4))) */
1394
1395 set_float_rounding_mode(user_rnd_mode, status);
1396 set_floatx80_rounding_precision(user_rnd_prec, status);
1397
1398 a = floatx80_div(fp0, fp1, status);
1399
1400 float_raise(float_flag_inexact, status);
1401
1402 return a;
1403 }
1404 }
1405 }
1406
1407 /*
1408 * Sine
1409 */
1410
1411 floatx80 floatx80_sin(floatx80 a, float_status *status)
1412 {
1413 bool aSign, xSign;
1414 int32_t aExp, xExp;
1415 uint64_t aSig, xSig;
1416
1417 FloatRoundMode user_rnd_mode;
1418 FloatX80RoundPrec user_rnd_prec;
1419
1420 int32_t compact, l, n, j;
1421 floatx80 fp0, fp1, fp2, fp3, fp4, fp5, x, invtwopi, twopi1, twopi2;
1422 float32 posneg1, twoto63;
1423 bool endflag;
1424
1425 aSig = extractFloatx80Frac(a);
1426 aExp = extractFloatx80Exp(a);
1427 aSign = extractFloatx80Sign(a);
1428
1429 if (aExp == 0x7FFF) {
1430 if ((uint64_t) (aSig << 1)) {
1431 return propagateFloatx80NaNOneArg(a, status);
1432 }
1433 float_raise(float_flag_invalid, status);
1434 return floatx80_default_nan(status);
1435 }
1436
1437 if (aExp == 0 && aSig == 0) {
1438 return packFloatx80(aSign, 0, 0);
1439 }
1440
1441 user_rnd_mode = get_float_rounding_mode(status);
1442 user_rnd_prec = get_floatx80_rounding_precision(status);
1443 set_float_rounding_mode(float_round_nearest_even, status);
1444 set_floatx80_rounding_precision(floatx80_precision_x, status);
1445
1446 compact = floatx80_make_compact(aExp, aSig);
1447
1448 fp0 = a;
1449
1450 if (compact < 0x3FD78000 || compact > 0x4004BC7E) {
1451 /* 2^(-40) > |X| > 15 PI */
1452 if (compact > 0x3FFF8000) { /* |X| >= 15 PI */
1453 /* REDUCEX */
1454 fp1 = packFloatx80(0, 0, 0);
1455 if (compact == 0x7FFEFFFF) {
1456 twopi1 = packFloatx80(aSign ^ 1, 0x7FFE,
1457 UINT64_C(0xC90FDAA200000000));
1458 twopi2 = packFloatx80(aSign ^ 1, 0x7FDC,
1459 UINT64_C(0x85A308D300000000));
1460 fp0 = floatx80_add(fp0, twopi1, status);
1461 fp1 = fp0;
1462 fp0 = floatx80_add(fp0, twopi2, status);
1463 fp1 = floatx80_sub(fp1, fp0, status);
1464 fp1 = floatx80_add(fp1, twopi2, status);
1465 }
1466 loop:
1467 xSign = extractFloatx80Sign(fp0);
1468 xExp = extractFloatx80Exp(fp0);
1469 xExp -= 0x3FFF;
1470 if (xExp <= 28) {
1471 l = 0;
1472 endflag = true;
1473 } else {
1474 l = xExp - 27;
1475 endflag = false;
1476 }
1477 invtwopi = packFloatx80(0, 0x3FFE - l,
1478 UINT64_C(0xA2F9836E4E44152A)); /* INVTWOPI */
1479 twopi1 = packFloatx80(0, 0x3FFF + l, UINT64_C(0xC90FDAA200000000));
1480 twopi2 = packFloatx80(0, 0x3FDD + l, UINT64_C(0x85A308D300000000));
1481
1482 /* SIGN(INARG)*2^63 IN SGL */
1483 twoto63 = packFloat32(xSign, 0xBE, 0);
1484
1485 fp2 = floatx80_mul(fp0, invtwopi, status);
1486 fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63, status),
1487 status); /* THE FRACT PART OF FP2 IS ROUNDED */
1488 fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63, status),
1489 status); /* FP2 is N */
1490 fp4 = floatx80_mul(twopi1, fp2, status); /* W = N*P1 */
1491 fp5 = floatx80_mul(twopi2, fp2, status); /* w = N*P2 */
1492 fp3 = floatx80_add(fp4, fp5, status); /* FP3 is P */
1493 fp4 = floatx80_sub(fp4, fp3, status); /* W-P */
1494 fp0 = floatx80_sub(fp0, fp3, status); /* FP0 is A := R - P */
1495 fp4 = floatx80_add(fp4, fp5, status); /* FP4 is p = (W-P)+w */
1496 fp3 = fp0; /* FP3 is A */
1497 fp1 = floatx80_sub(fp1, fp4, status); /* FP1 is a := r - p */
1498 fp0 = floatx80_add(fp0, fp1, status); /* FP0 is R := A+a */
1499
1500 if (endflag) {
1501 n = floatx80_to_int32(fp2, status);
1502 goto sincont;
1503 }
1504 fp3 = floatx80_sub(fp3, fp0, status); /* A-R */
1505 fp1 = floatx80_add(fp1, fp3, status); /* FP1 is r := (A-R)+a */
1506 goto loop;
1507 } else {
1508 /* SINSM */
1509 fp0 = float32_to_floatx80(make_float32(0x3F800000),
1510 status); /* 1 */
1511
1512 set_float_rounding_mode(user_rnd_mode, status);
1513 set_floatx80_rounding_precision(user_rnd_prec, status);
1514
1515 /* SINTINY */
1516 a = floatx80_move(a, status);
1517 float_raise(float_flag_inexact, status);
1518
1519 return a;
1520 }
1521 } else {
1522 fp1 = floatx80_mul(fp0, float64_to_floatx80(
1523 make_float64(0x3FE45F306DC9C883), status),
1524 status); /* X*2/PI */
1525
1526 n = floatx80_to_int32(fp1, status);
1527 j = 32 + n;
1528
1529 fp0 = floatx80_sub(fp0, pi_tbl[j], status); /* X-Y1 */
1530 fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j], status),
1531 status); /* FP0 IS R = (X-Y1)-Y2 */
1532
1533 sincont:
1534 if (n & 1) {
1535 /* COSPOLY */
1536 fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */
1537 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */
1538 fp2 = float64_to_floatx80(make_float64(0x3D2AC4D0D6011EE3),
1539 status); /* B8 */
1540 fp3 = float64_to_floatx80(make_float64(0xBDA9396F9F45AC19),
1541 status); /* B7 */
1542
1543 xSign = extractFloatx80Sign(fp0); /* X IS S */
1544 xExp = extractFloatx80Exp(fp0);
1545 xSig = extractFloatx80Frac(fp0);
1546
1547 if ((n >> 1) & 1) {
1548 xSign ^= 1;
1549 posneg1 = make_float32(0xBF800000); /* -1 */
1550 } else {
1551 xSign ^= 0;
1552 posneg1 = make_float32(0x3F800000); /* 1 */
1553 } /* X IS NOW R'= SGN*R */
1554
1555 fp2 = floatx80_mul(fp2, fp1, status); /* TB8 */
1556 fp3 = floatx80_mul(fp3, fp1, status); /* TB7 */
1557 fp2 = floatx80_add(fp2, float64_to_floatx80(
1558 make_float64(0x3E21EED90612C972), status),
1559 status); /* B6+TB8 */
1560 fp3 = floatx80_add(fp3, float64_to_floatx80(
1561 make_float64(0xBE927E4FB79D9FCF), status),
1562 status); /* B5+TB7 */
1563 fp2 = floatx80_mul(fp2, fp1, status); /* T(B6+TB8) */
1564 fp3 = floatx80_mul(fp3, fp1, status); /* T(B5+TB7) */
1565 fp2 = floatx80_add(fp2, float64_to_floatx80(
1566 make_float64(0x3EFA01A01A01D423), status),
1567 status); /* B4+T(B6+TB8) */
1568 fp4 = packFloatx80(1, 0x3FF5, UINT64_C(0xB60B60B60B61D438));
1569 fp3 = floatx80_add(fp3, fp4, status); /* B3+T(B5+TB7) */
1570 fp2 = floatx80_mul(fp2, fp1, status); /* T(B4+T(B6+TB8)) */
1571 fp1 = floatx80_mul(fp1, fp3, status); /* T(B3+T(B5+TB7)) */
1572 fp4 = packFloatx80(0, 0x3FFA, UINT64_C(0xAAAAAAAAAAAAAB5E));
1573 fp2 = floatx80_add(fp2, fp4, status); /* B2+T(B4+T(B6+TB8)) */
1574 fp1 = floatx80_add(fp1, float32_to_floatx80(
1575 make_float32(0xBF000000), status),
1576 status); /* B1+T(B3+T(B5+TB7)) */
1577 fp0 = floatx80_mul(fp0, fp2, status); /* S(B2+T(B4+T(B6+TB8))) */
1578 fp0 = floatx80_add(fp0, fp1, status); /* [B1+T(B3+T(B5+TB7))]+
1579 * [S(B2+T(B4+T(B6+TB8)))]
1580 */
1581
1582 x = packFloatx80(xSign, xExp, xSig);
1583 fp0 = floatx80_mul(fp0, x, status);
1584
1585 set_float_rounding_mode(user_rnd_mode, status);
1586 set_floatx80_rounding_precision(user_rnd_prec, status);
1587
1588 a = floatx80_add(fp0, float32_to_floatx80(posneg1, status), status);
1589
1590 float_raise(float_flag_inexact, status);
1591
1592 return a;
1593 } else {
1594 /* SINPOLY */
1595 xSign = extractFloatx80Sign(fp0); /* X IS R */
1596 xExp = extractFloatx80Exp(fp0);
1597 xSig = extractFloatx80Frac(fp0);
1598
1599 xSign ^= (n >> 1) & 1; /* X IS NOW R'= SGN*R */
1600
1601 fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */
1602 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */
1603 fp3 = float64_to_floatx80(make_float64(0xBD6AAA77CCC994F5),
1604 status); /* A7 */
1605 fp2 = float64_to_floatx80(make_float64(0x3DE612097AAE8DA1),
1606 status); /* A6 */
1607 fp3 = floatx80_mul(fp3, fp1, status); /* T*A7 */
1608 fp2 = floatx80_mul(fp2, fp1, status); /* T*A6 */
1609 fp3 = floatx80_add(fp3, float64_to_floatx80(
1610 make_float64(0xBE5AE6452A118AE4), status),
1611 status); /* A5+T*A7 */
1612 fp2 = floatx80_add(fp2, float64_to_floatx80(
1613 make_float64(0x3EC71DE3A5341531), status),
1614 status); /* A4+T*A6 */
1615 fp3 = floatx80_mul(fp3, fp1, status); /* T(A5+TA7) */
1616 fp2 = floatx80_mul(fp2, fp1, status); /* T(A4+TA6) */
1617 fp3 = floatx80_add(fp3, float64_to_floatx80(
1618 make_float64(0xBF2A01A01A018B59), status),
1619 status); /* A3+T(A5+TA7) */
1620 fp4 = packFloatx80(0, 0x3FF8, UINT64_C(0x88888888888859AF));
1621 fp2 = floatx80_add(fp2, fp4, status); /* A2+T(A4+TA6) */
1622 fp1 = floatx80_mul(fp1, fp3, status); /* T(A3+T(A5+TA7)) */
1623 fp2 = floatx80_mul(fp2, fp0, status); /* S(A2+T(A4+TA6)) */
1624 fp4 = packFloatx80(1, 0x3FFC, UINT64_C(0xAAAAAAAAAAAAAA99));
1625 fp1 = floatx80_add(fp1, fp4, status); /* A1+T(A3+T(A5+TA7)) */
1626 fp1 = floatx80_add(fp1, fp2,
1627 status); /* [A1+T(A3+T(A5+TA7))]+
1628 * [S(A2+T(A4+TA6))]
1629 */
1630
1631 x = packFloatx80(xSign, xExp, xSig);
1632 fp0 = floatx80_mul(fp0, x, status); /* R'*S */
1633 fp0 = floatx80_mul(fp0, fp1, status); /* SIN(R')-R' */
1634
1635 set_float_rounding_mode(user_rnd_mode, status);
1636 set_floatx80_rounding_precision(user_rnd_prec, status);
1637
1638 a = floatx80_add(fp0, x, status);
1639
1640 float_raise(float_flag_inexact, status);
1641
1642 return a;
1643 }
1644 }
1645 }
1646
1647 /*
1648 * Cosine
1649 */
1650
1651 floatx80 floatx80_cos(floatx80 a, float_status *status)
1652 {
1653 bool aSign, xSign;
1654 int32_t aExp, xExp;
1655 uint64_t aSig, xSig;
1656
1657 FloatRoundMode user_rnd_mode;
1658 FloatX80RoundPrec user_rnd_prec;
1659
1660 int32_t compact, l, n, j;
1661 floatx80 fp0, fp1, fp2, fp3, fp4, fp5, x, invtwopi, twopi1, twopi2;
1662 float32 posneg1, twoto63;
1663 bool endflag;
1664
1665 aSig = extractFloatx80Frac(a);
1666 aExp = extractFloatx80Exp(a);
1667 aSign = extractFloatx80Sign(a);
1668
1669 if (aExp == 0x7FFF) {
1670 if ((uint64_t) (aSig << 1)) {
1671 return propagateFloatx80NaNOneArg(a, status);
1672 }
1673 float_raise(float_flag_invalid, status);
1674 return floatx80_default_nan(status);
1675 }
1676
1677 if (aExp == 0 && aSig == 0) {
1678 return packFloatx80(0, one_exp, one_sig);
1679 }
1680
1681 user_rnd_mode = get_float_rounding_mode(status);
1682 user_rnd_prec = get_floatx80_rounding_precision(status);
1683 set_float_rounding_mode(float_round_nearest_even, status);
1684 set_floatx80_rounding_precision(floatx80_precision_x, status);
1685
1686 compact = floatx80_make_compact(aExp, aSig);
1687
1688 fp0 = a;
1689
1690 if (compact < 0x3FD78000 || compact > 0x4004BC7E) {
1691 /* 2^(-40) > |X| > 15 PI */
1692 if (compact > 0x3FFF8000) { /* |X| >= 15 PI */
1693 /* REDUCEX */
1694 fp1 = packFloatx80(0, 0, 0);
1695 if (compact == 0x7FFEFFFF) {
1696 twopi1 = packFloatx80(aSign ^ 1, 0x7FFE,
1697 UINT64_C(0xC90FDAA200000000));
1698 twopi2 = packFloatx80(aSign ^ 1, 0x7FDC,
1699 UINT64_C(0x85A308D300000000));
1700 fp0 = floatx80_add(fp0, twopi1, status);
1701 fp1 = fp0;
1702 fp0 = floatx80_add(fp0, twopi2, status);
1703 fp1 = floatx80_sub(fp1, fp0, status);
1704 fp1 = floatx80_add(fp1, twopi2, status);
1705 }
1706 loop:
1707 xSign = extractFloatx80Sign(fp0);
1708 xExp = extractFloatx80Exp(fp0);
1709 xExp -= 0x3FFF;
1710 if (xExp <= 28) {
1711 l = 0;
1712 endflag = true;
1713 } else {
1714 l = xExp - 27;
1715 endflag = false;
1716 }
1717 invtwopi = packFloatx80(0, 0x3FFE - l,
1718 UINT64_C(0xA2F9836E4E44152A)); /* INVTWOPI */
1719 twopi1 = packFloatx80(0, 0x3FFF + l, UINT64_C(0xC90FDAA200000000));
1720 twopi2 = packFloatx80(0, 0x3FDD + l, UINT64_C(0x85A308D300000000));
1721
1722 /* SIGN(INARG)*2^63 IN SGL */
1723 twoto63 = packFloat32(xSign, 0xBE, 0);
1724
1725 fp2 = floatx80_mul(fp0, invtwopi, status);
1726 fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63, status),
1727 status); /* THE FRACT PART OF FP2 IS ROUNDED */
1728 fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63, status),
1729 status); /* FP2 is N */
1730 fp4 = floatx80_mul(twopi1, fp2, status); /* W = N*P1 */
1731 fp5 = floatx80_mul(twopi2, fp2, status); /* w = N*P2 */
1732 fp3 = floatx80_add(fp4, fp5, status); /* FP3 is P */
1733 fp4 = floatx80_sub(fp4, fp3, status); /* W-P */
1734 fp0 = floatx80_sub(fp0, fp3, status); /* FP0 is A := R - P */
1735 fp4 = floatx80_add(fp4, fp5, status); /* FP4 is p = (W-P)+w */
1736 fp3 = fp0; /* FP3 is A */
1737 fp1 = floatx80_sub(fp1, fp4, status); /* FP1 is a := r - p */
1738 fp0 = floatx80_add(fp0, fp1, status); /* FP0 is R := A+a */
1739
1740 if (endflag) {
1741 n = floatx80_to_int32(fp2, status);
1742 goto sincont;
1743 }
1744 fp3 = floatx80_sub(fp3, fp0, status); /* A-R */
1745 fp1 = floatx80_add(fp1, fp3, status); /* FP1 is r := (A-R)+a */
1746 goto loop;
1747 } else {
1748 /* SINSM */
1749 fp0 = float32_to_floatx80(make_float32(0x3F800000), status); /* 1 */
1750
1751 set_float_rounding_mode(user_rnd_mode, status);
1752 set_floatx80_rounding_precision(user_rnd_prec, status);
1753
1754 /* COSTINY */
1755 a = floatx80_sub(fp0, float32_to_floatx80(
1756 make_float32(0x00800000), status),
1757 status);
1758 float_raise(float_flag_inexact, status);
1759
1760 return a;
1761 }
1762 } else {
1763 fp1 = floatx80_mul(fp0, float64_to_floatx80(
1764 make_float64(0x3FE45F306DC9C883), status),
1765 status); /* X*2/PI */
1766
1767 n = floatx80_to_int32(fp1, status);
1768 j = 32 + n;
1769
1770 fp0 = floatx80_sub(fp0, pi_tbl[j], status); /* X-Y1 */
1771 fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j], status),
1772 status); /* FP0 IS R = (X-Y1)-Y2 */
1773
1774 sincont:
1775 if ((n + 1) & 1) {
1776 /* COSPOLY */
1777 fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */
1778 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */
1779 fp2 = float64_to_floatx80(make_float64(0x3D2AC4D0D6011EE3),
1780 status); /* B8 */
1781 fp3 = float64_to_floatx80(make_float64(0xBDA9396F9F45AC19),
1782 status); /* B7 */
1783
1784 xSign = extractFloatx80Sign(fp0); /* X IS S */
1785 xExp = extractFloatx80Exp(fp0);
1786 xSig = extractFloatx80Frac(fp0);
1787
1788 if (((n + 1) >> 1) & 1) {
1789 xSign ^= 1;
1790 posneg1 = make_float32(0xBF800000); /* -1 */
1791 } else {
1792 xSign ^= 0;
1793 posneg1 = make_float32(0x3F800000); /* 1 */
1794 } /* X IS NOW R'= SGN*R */
1795
1796 fp2 = floatx80_mul(fp2, fp1, status); /* TB8 */
1797 fp3 = floatx80_mul(fp3, fp1, status); /* TB7 */
1798 fp2 = floatx80_add(fp2, float64_to_floatx80(
1799 make_float64(0x3E21EED90612C972), status),
1800 status); /* B6+TB8 */
1801 fp3 = floatx80_add(fp3, float64_to_floatx80(
1802 make_float64(0xBE927E4FB79D9FCF), status),
1803 status); /* B5+TB7 */
1804 fp2 = floatx80_mul(fp2, fp1, status); /* T(B6+TB8) */
1805 fp3 = floatx80_mul(fp3, fp1, status); /* T(B5+TB7) */
1806 fp2 = floatx80_add(fp2, float64_to_floatx80(
1807 make_float64(0x3EFA01A01A01D423), status),
1808 status); /* B4+T(B6+TB8) */
1809 fp4 = packFloatx80(1, 0x3FF5, UINT64_C(0xB60B60B60B61D438));
1810 fp3 = floatx80_add(fp3, fp4, status); /* B3+T(B5+TB7) */
1811 fp2 = floatx80_mul(fp2, fp1, status); /* T(B4+T(B6+TB8)) */
1812 fp1 = floatx80_mul(fp1, fp3, status); /* T(B3+T(B5+TB7)) */
1813 fp4 = packFloatx80(0, 0x3FFA, UINT64_C(0xAAAAAAAAAAAAAB5E));
1814 fp2 = floatx80_add(fp2, fp4, status); /* B2+T(B4+T(B6+TB8)) */
1815 fp1 = floatx80_add(fp1, float32_to_floatx80(
1816 make_float32(0xBF000000), status),
1817 status); /* B1+T(B3+T(B5+TB7)) */
1818 fp0 = floatx80_mul(fp0, fp2, status); /* S(B2+T(B4+T(B6+TB8))) */
1819 fp0 = floatx80_add(fp0, fp1, status);
1820 /* [B1+T(B3+T(B5+TB7))]+[S(B2+T(B4+T(B6+TB8)))] */
1821
1822 x = packFloatx80(xSign, xExp, xSig);
1823 fp0 = floatx80_mul(fp0, x, status);
1824
1825 set_float_rounding_mode(user_rnd_mode, status);
1826 set_floatx80_rounding_precision(user_rnd_prec, status);
1827
1828 a = floatx80_add(fp0, float32_to_floatx80(posneg1, status), status);
1829
1830 float_raise(float_flag_inexact, status);
1831
1832 return a;
1833 } else {
1834 /* SINPOLY */
1835 xSign = extractFloatx80Sign(fp0); /* X IS R */
1836 xExp = extractFloatx80Exp(fp0);
1837 xSig = extractFloatx80Frac(fp0);
1838
1839 xSign ^= ((n + 1) >> 1) & 1; /* X IS NOW R'= SGN*R */
1840
1841 fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */
1842 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */
1843 fp3 = float64_to_floatx80(make_float64(0xBD6AAA77CCC994F5),
1844 status); /* A7 */
1845 fp2 = float64_to_floatx80(make_float64(0x3DE612097AAE8DA1),
1846 status); /* A6 */
1847 fp3 = floatx80_mul(fp3, fp1, status); /* T*A7 */
1848 fp2 = floatx80_mul(fp2, fp1, status); /* T*A6 */
1849 fp3 = floatx80_add(fp3, float64_to_floatx80(
1850 make_float64(0xBE5AE6452A118AE4), status),
1851 status); /* A5+T*A7 */
1852 fp2 = floatx80_add(fp2, float64_to_floatx80(
1853 make_float64(0x3EC71DE3A5341531), status),
1854 status); /* A4+T*A6 */
1855 fp3 = floatx80_mul(fp3, fp1, status); /* T(A5+TA7) */
1856 fp2 = floatx80_mul(fp2, fp1, status); /* T(A4+TA6) */
1857 fp3 = floatx80_add(fp3, float64_to_floatx80(
1858 make_float64(0xBF2A01A01A018B59), status),
1859 status); /* A3+T(A5+TA7) */
1860 fp4 = packFloatx80(0, 0x3FF8, UINT64_C(0x88888888888859AF));
1861 fp2 = floatx80_add(fp2, fp4, status); /* A2+T(A4+TA6) */
1862 fp1 = floatx80_mul(fp1, fp3, status); /* T(A3+T(A5+TA7)) */
1863 fp2 = floatx80_mul(fp2, fp0, status); /* S(A2+T(A4+TA6)) */
1864 fp4 = packFloatx80(1, 0x3FFC, UINT64_C(0xAAAAAAAAAAAAAA99));
1865 fp1 = floatx80_add(fp1, fp4, status); /* A1+T(A3+T(A5+TA7)) */
1866 fp1 = floatx80_add(fp1, fp2, status);
1867 /* [A1+T(A3+T(A5+TA7))]+[S(A2+T(A4+TA6))] */
1868
1869 x = packFloatx80(xSign, xExp, xSig);
1870 fp0 = floatx80_mul(fp0, x, status); /* R'*S */
1871 fp0 = floatx80_mul(fp0, fp1, status); /* SIN(R')-R' */
1872
1873 set_float_rounding_mode(user_rnd_mode, status);
1874 set_floatx80_rounding_precision(user_rnd_prec, status);
1875
1876 a = floatx80_add(fp0, x, status);
1877
1878 float_raise(float_flag_inexact, status);
1879
1880 return a;
1881 }
1882 }
1883 }
1884
1885 /*
1886 * Arc tangent
1887 */
1888
1889 floatx80 floatx80_atan(floatx80 a, float_status *status)
1890 {
1891 bool aSign;
1892 int32_t aExp;
1893 uint64_t aSig;
1894
1895 FloatRoundMode user_rnd_mode;
1896 FloatX80RoundPrec user_rnd_prec;
1897
1898 int32_t compact, tbl_index;
1899 floatx80 fp0, fp1, fp2, fp3, xsave;
1900
1901 aSig = extractFloatx80Frac(a);
1902 aExp = extractFloatx80Exp(a);
1903 aSign = extractFloatx80Sign(a);
1904
1905 if (aExp == 0x7FFF) {
1906 if ((uint64_t) (aSig << 1)) {
1907 return propagateFloatx80NaNOneArg(a, status);
1908 }
1909 a = packFloatx80(aSign, piby2_exp, pi_sig);
1910 float_raise(float_flag_inexact, status);
1911 return floatx80_move(a, status);
1912 }
1913
1914 if (aExp == 0 && aSig == 0) {
1915 return packFloatx80(aSign, 0, 0);
1916 }
1917
1918 compact = floatx80_make_compact(aExp, aSig);
1919
1920 user_rnd_mode = get_float_rounding_mode(status);
1921 user_rnd_prec = get_floatx80_rounding_precision(status);
1922 set_float_rounding_mode(float_round_nearest_even, status);
1923 set_floatx80_rounding_precision(floatx80_precision_x, status);
1924
1925 if (compact < 0x3FFB8000 || compact > 0x4002FFFF) {
1926 /* |X| >= 16 or |X| < 1/16 */
1927 if (compact > 0x3FFF8000) { /* |X| >= 16 */
1928 if (compact > 0x40638000) { /* |X| > 2^(100) */
1929 fp0 = packFloatx80(aSign, piby2_exp, pi_sig);
1930 fp1 = packFloatx80(aSign, 0x0001, one_sig);
1931
1932 set_float_rounding_mode(user_rnd_mode, status);
1933 set_floatx80_rounding_precision(user_rnd_prec, status);
1934
1935 a = floatx80_sub(fp0, fp1, status);
1936
1937 float_raise(float_flag_inexact, status);
1938
1939 return a;
1940 } else {
1941 fp0 = a;
1942 fp1 = packFloatx80(1, one_exp, one_sig); /* -1 */
1943 fp1 = floatx80_div(fp1, fp0, status); /* X' = -1/X */
1944 xsave = fp1;
1945 fp0 = floatx80_mul(fp1, fp1, status); /* Y = X'*X' */
1946 fp1 = floatx80_mul(fp0, fp0, status); /* Z = Y*Y */
1947 fp3 = float64_to_floatx80(make_float64(0xBFB70BF398539E6A),
1948 status); /* C5 */
1949 fp2 = float64_to_floatx80(make_float64(0x3FBC7187962D1D7D),
1950 status); /* C4 */
1951 fp3 = floatx80_mul(fp3, fp1, status); /* Z*C5 */
1952 fp2 = floatx80_mul(fp2, fp1, status); /* Z*C4 */
1953 fp3 = floatx80_add(fp3, float64_to_floatx80(
1954 make_float64(0xBFC24924827107B8), status),
1955 status); /* C3+Z*C5 */
1956 fp2 = floatx80_add(fp2, float64_to_floatx80(
1957 make_float64(0x3FC999999996263E), status),
1958 status); /* C2+Z*C4 */
1959 fp1 = floatx80_mul(fp1, fp3, status); /* Z*(C3+Z*C5) */
1960 fp2 = floatx80_mul(fp2, fp0, status); /* Y*(C2+Z*C4) */
1961 fp1 = floatx80_add(fp1, float64_to_floatx80(
1962 make_float64(0xBFD5555555555536), status),
1963 status); /* C1+Z*(C3+Z*C5) */
1964 fp0 = floatx80_mul(fp0, xsave, status); /* X'*Y */
1965 /* [Y*(C2+Z*C4)]+[C1+Z*(C3+Z*C5)] */
1966 fp1 = floatx80_add(fp1, fp2, status);
1967 /* X'*Y*([B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))]) ?? */
1968 fp0 = floatx80_mul(fp0, fp1, status);
1969 fp0 = floatx80_add(fp0, xsave, status);
1970 fp1 = packFloatx80(aSign, piby2_exp, pi_sig);
1971
1972 set_float_rounding_mode(user_rnd_mode, status);
1973 set_floatx80_rounding_precision(user_rnd_prec, status);
1974
1975 a = floatx80_add(fp0, fp1, status);
1976
1977 float_raise(float_flag_inexact, status);
1978
1979 return a;
1980 }
1981 } else { /* |X| < 1/16 */
1982 if (compact < 0x3FD78000) { /* |X| < 2^(-40) */
1983 set_float_rounding_mode(user_rnd_mode, status);
1984 set_floatx80_rounding_precision(user_rnd_prec, status);
1985
1986 a = floatx80_move(a, status);
1987
1988 float_raise(float_flag_inexact, status);
1989
1990 return a;
1991 } else {
1992 fp0 = a;
1993 xsave = a;
1994 fp0 = floatx80_mul(fp0, fp0, status); /* Y = X*X */
1995 fp1 = floatx80_mul(fp0, fp0, status); /* Z = Y*Y */
1996 fp2 = float64_to_floatx80(make_float64(0x3FB344447F876989),
1997 status); /* B6 */
1998 fp3 = float64_to_floatx80(make_float64(0xBFB744EE7FAF45DB),
1999 status); /* B5 */
2000 fp2 = floatx80_mul(fp2, fp1, status); /* Z*B6 */
2001 fp3 = floatx80_mul(fp3, fp1, status); /* Z*B5 */
2002 fp2 = floatx80_add(fp2, float64_to_floatx80(
2003 make_float64(0x3FBC71C646940220), status),
2004 status); /* B4+Z*B6 */
2005 fp3 = floatx80_add(fp3, float64_to_floatx80(
2006 make_float64(0xBFC24924921872F9),
2007 status), status); /* B3+Z*B5 */
2008 fp2 = floatx80_mul(fp2, fp1, status); /* Z*(B4+Z*B6) */
2009 fp1 = floatx80_mul(fp1, fp3, status); /* Z*(B3+Z*B5) */
2010 fp2 = floatx80_add(fp2, float64_to_floatx80(
2011 make_float64(0x3FC9999999998FA9), status),
2012 status); /* B2+Z*(B4+Z*B6) */
2013 fp1 = floatx80_add(fp1, float64_to_floatx80(
2014 make_float64(0xBFD5555555555555), status),
2015 status); /* B1+Z*(B3+Z*B5) */
2016 fp2 = floatx80_mul(fp2, fp0, status); /* Y*(B2+Z*(B4+Z*B6)) */
2017 fp0 = floatx80_mul(fp0, xsave, status); /* X*Y */
2018 /* [B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))] */
2019 fp1 = floatx80_add(fp1, fp2, status);
2020 /* X*Y*([B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))]) */
2021 fp0 = floatx80_mul(fp0, fp1, status);
2022
2023 set_float_rounding_mode(user_rnd_mode, status);
2024 set_floatx80_rounding_precision(user_rnd_prec, status);
2025
2026 a = floatx80_add(fp0, xsave, status);
2027
2028 float_raise(float_flag_inexact, status);
2029
2030 return a;
2031 }
2032 }
2033 } else {
2034 aSig &= UINT64_C(0xF800000000000000);
2035 aSig |= UINT64_C(0x0400000000000000);
2036 xsave = packFloatx80(aSign, aExp, aSig); /* F */
2037 fp0 = a;
2038 fp1 = a; /* X */
2039 fp2 = packFloatx80(0, one_exp, one_sig); /* 1 */
2040 fp1 = floatx80_mul(fp1, xsave, status); /* X*F */
2041 fp0 = floatx80_sub(fp0, xsave, status); /* X-F */
2042 fp1 = floatx80_add(fp1, fp2, status); /* 1 + X*F */
2043 fp0 = floatx80_div(fp0, fp1, status); /* U = (X-F)/(1+X*F) */
2044
2045 tbl_index = compact;
2046
2047 tbl_index &= 0x7FFF0000;
2048 tbl_index -= 0x3FFB0000;
2049 tbl_index >>= 1;
2050 tbl_index += compact & 0x00007800;
2051 tbl_index >>= 11;
2052
2053 fp3 = atan_tbl[tbl_index];
2054
2055 fp3.high |= aSign ? 0x8000 : 0; /* ATAN(F) */
2056
2057 fp1 = floatx80_mul(fp0, fp0, status); /* V = U*U */
2058 fp2 = float64_to_floatx80(make_float64(0xBFF6687E314987D8),
2059 status); /* A3 */
2060 fp2 = floatx80_add(fp2, fp1, status); /* A3+V */
2061 fp2 = floatx80_mul(fp2, fp1, status); /* V*(A3+V) */
2062 fp1 = floatx80_mul(fp1, fp0, status); /* U*V */
2063 fp2 = floatx80_add(fp2, float64_to_floatx80(
2064 make_float64(0x4002AC6934A26DB3), status),
2065 status); /* A2+V*(A3+V) */
2066 fp1 = floatx80_mul(fp1, float64_to_floatx80(
2067 make_float64(0xBFC2476F4E1DA28E), status),
2068 status); /* A1+U*V */
2069 fp1 = floatx80_mul(fp1, fp2, status); /* A1*U*V*(A2+V*(A3+V)) */
2070 fp0 = floatx80_add(fp0, fp1, status); /* ATAN(U) */
2071
2072 set_float_rounding_mode(user_rnd_mode, status);
2073 set_floatx80_rounding_precision(user_rnd_prec, status);
2074
2075 a = floatx80_add(fp0, fp3, status); /* ATAN(X) */
2076
2077 float_raise(float_flag_inexact, status);
2078
2079 return a;
2080 }
2081 }
2082
2083 /*
2084 * Arc sine
2085 */
2086
2087 floatx80 floatx80_asin(floatx80 a, float_status *status)
2088 {
2089 bool aSign;
2090 int32_t aExp;
2091 uint64_t aSig;
2092
2093 FloatRoundMode user_rnd_mode;
2094 FloatX80RoundPrec user_rnd_prec;
2095
2096 int32_t compact;
2097 floatx80 fp0, fp1, fp2, one;
2098
2099 aSig = extractFloatx80Frac(a);
2100 aExp = extractFloatx80Exp(a);
2101 aSign = extractFloatx80Sign(a);
2102
2103 if (aExp == 0x7FFF && (uint64_t) (aSig << 1)) {
2104 return propagateFloatx80NaNOneArg(a, status);
2105 }
2106
2107 if (aExp == 0 && aSig == 0) {
2108 return packFloatx80(aSign, 0, 0);
2109 }
2110
2111 compact = floatx80_make_compact(aExp, aSig);
2112
2113 if (compact >= 0x3FFF8000) { /* |X| >= 1 */
2114 if (aExp == one_exp && aSig == one_sig) { /* |X| == 1 */
2115 float_raise(float_flag_inexact, status);
2116 a = packFloatx80(aSign, piby2_exp, pi_sig);
2117 return floatx80_move(a, status);
2118 } else { /* |X| > 1 */
2119 float_raise(float_flag_invalid, status);
2120 return floatx80_default_nan(status);
2121 }
2122
2123 } /* |X| < 1 */
2124
2125 user_rnd_mode = get_float_rounding_mode(status);
2126 user_rnd_prec = get_floatx80_rounding_precision(status);
2127 set_float_rounding_mode(float_round_nearest_even, status);
2128 set_floatx80_rounding_precision(floatx80_precision_x, status);
2129
2130 one = packFloatx80(0, one_exp, one_sig);
2131 fp0 = a;
2132
2133 fp1 = floatx80_sub(one, fp0, status); /* 1 - X */
2134 fp2 = floatx80_add(one, fp0, status); /* 1 + X */
2135 fp1 = floatx80_mul(fp2, fp1, status); /* (1+X)*(1-X) */
2136 fp1 = floatx80_sqrt(fp1, status); /* SQRT((1+X)*(1-X)) */
2137 fp0 = floatx80_div(fp0, fp1, status); /* X/SQRT((1+X)*(1-X)) */
2138
2139 set_float_rounding_mode(user_rnd_mode, status);
2140 set_floatx80_rounding_precision(user_rnd_prec, status);
2141
2142 a = floatx80_atan(fp0, status); /* ATAN(X/SQRT((1+X)*(1-X))) */
2143
2144 float_raise(float_flag_inexact, status);
2145
2146 return a;
2147 }
2148
2149 /*
2150 * Arc cosine
2151 */
2152
2153 floatx80 floatx80_acos(floatx80 a, float_status *status)
2154 {
2155 bool aSign;
2156 int32_t aExp;
2157 uint64_t aSig;
2158
2159 FloatRoundMode user_rnd_mode;
2160 FloatX80RoundPrec user_rnd_prec;
2161
2162 int32_t compact;
2163 floatx80 fp0, fp1, one;
2164
2165 aSig = extractFloatx80Frac(a);
2166 aExp = extractFloatx80Exp(a);
2167 aSign = extractFloatx80Sign(a);
2168
2169 if (aExp == 0x7FFF && (uint64_t) (aSig << 1)) {
2170 return propagateFloatx80NaNOneArg(a, status);
2171 }
2172 if (aExp == 0 && aSig == 0) {
2173 float_raise(float_flag_inexact, status);
2174 return roundAndPackFloatx80(get_floatx80_rounding_precision(status),
2175 0,
2176 piby2_exp, pi_sig, 0, status);
2177 }
2178
2179 compact = floatx80_make_compact(aExp, aSig);
2180
2181 if (compact >= 0x3FFF8000) { /* |X| >= 1 */
2182 if (aExp == one_exp && aSig == one_sig) { /* |X| == 1 */
2183 if (aSign) { /* X == -1 */
2184 a = packFloatx80(0, pi_exp, pi_sig);
2185 float_raise(float_flag_inexact, status);
2186 return floatx80_move(a, status);
2187 } else { /* X == +1 */
2188 return packFloatx80(0, 0, 0);
2189 }
2190 } else { /* |X| > 1 */
2191 float_raise(float_flag_invalid, status);
2192 return floatx80_default_nan(status);
2193 }
2194 } /* |X| < 1 */
2195
2196 user_rnd_mode = get_float_rounding_mode(status);
2197 user_rnd_prec = get_floatx80_rounding_precision(status);
2198 set_float_rounding_mode(float_round_nearest_even, status);
2199 set_floatx80_rounding_precision(floatx80_precision_x, status);
2200
2201 one = packFloatx80(0, one_exp, one_sig);
2202 fp0 = a;
2203
2204 fp1 = floatx80_add(one, fp0, status); /* 1 + X */
2205 fp0 = floatx80_sub(one, fp0, status); /* 1 - X */
2206 fp0 = floatx80_div(fp0, fp1, status); /* (1-X)/(1+X) */
2207 fp0 = floatx80_sqrt(fp0, status); /* SQRT((1-X)/(1+X)) */
2208 fp0 = floatx80_atan(fp0, status); /* ATAN(SQRT((1-X)/(1+X))) */
2209
2210 set_float_rounding_mode(user_rnd_mode, status);
2211 set_floatx80_rounding_precision(user_rnd_prec, status);
2212
2213 a = floatx80_add(fp0, fp0, status); /* 2 * ATAN(SQRT((1-X)/(1+X))) */
2214
2215 float_raise(float_flag_inexact, status);
2216
2217 return a;
2218 }
2219
2220 /*
2221 * Hyperbolic arc tangent
2222 */
2223
2224 floatx80 floatx80_atanh(floatx80 a, float_status *status)
2225 {
2226 bool aSign;
2227 int32_t aExp;
2228 uint64_t aSig;
2229
2230 FloatRoundMode user_rnd_mode;
2231 FloatX80RoundPrec user_rnd_prec;
2232
2233 int32_t compact;
2234 floatx80 fp0, fp1, fp2, one;
2235
2236 aSig = extractFloatx80Frac(a);
2237 aExp = extractFloatx80Exp(a);
2238 aSign = extractFloatx80Sign(a);
2239
2240 if (aExp == 0x7FFF && (uint64_t) (aSig << 1)) {
2241 return propagateFloatx80NaNOneArg(a, status);
2242 }
2243
2244 if (aExp == 0 && aSig == 0) {
2245 return packFloatx80(aSign, 0, 0);
2246 }
2247
2248 compact = floatx80_make_compact(aExp, aSig);
2249
2250 if (compact >= 0x3FFF8000) { /* |X| >= 1 */
2251 if (aExp == one_exp && aSig == one_sig) { /* |X| == 1 */
2252 float_raise(float_flag_divbyzero, status);
2253 return floatx80_default_inf(aSign, status);
2254 } else { /* |X| > 1 */
2255 float_raise(float_flag_invalid, status);
2256 return floatx80_default_nan(status);
2257 }
2258 } /* |X| < 1 */
2259
2260 user_rnd_mode = get_float_rounding_mode(status);
2261 user_rnd_prec = get_floatx80_rounding_precision(status);
2262 set_float_rounding_mode(float_round_nearest_even, status);
2263 set_floatx80_rounding_precision(floatx80_precision_x, status);
2264
2265 one = packFloatx80(0, one_exp, one_sig);
2266 fp2 = packFloatx80(aSign, 0x3FFE, one_sig); /* SIGN(X) * (1/2) */
2267 fp0 = packFloatx80(0, aExp, aSig); /* Y = |X| */
2268 fp1 = packFloatx80(1, aExp, aSig); /* -Y */
2269 fp0 = floatx80_add(fp0, fp0, status); /* 2Y */
2270 fp1 = floatx80_add(fp1, one, status); /* 1-Y */
2271 fp0 = floatx80_div(fp0, fp1, status); /* Z = 2Y/(1-Y) */
2272 fp0 = floatx80_lognp1(fp0, status); /* LOG1P(Z) */
2273
2274 set_float_rounding_mode(user_rnd_mode, status);
2275 set_floatx80_rounding_precision(user_rnd_prec, status);
2276
2277 a = floatx80_mul(fp0, fp2,
2278 status); /* ATANH(X) = SIGN(X) * (1/2) * LOG1P(Z) */
2279
2280 float_raise(float_flag_inexact, status);
2281
2282 return a;
2283 }
2284
2285 /*
2286 * e to x minus 1
2287 */
2288
2289 floatx80 floatx80_etoxm1(floatx80 a, float_status *status)
2290 {
2291 bool aSign;
2292 int32_t aExp;
2293 uint64_t aSig;
2294
2295 FloatRoundMode user_rnd_mode;
2296 FloatX80RoundPrec user_rnd_prec;
2297
2298 int32_t compact, n, j, m, m1;
2299 floatx80 fp0, fp1, fp2, fp3, l2, sc, onebysc;
2300
2301 aSig = extractFloatx80Frac(a);
2302 aExp = extractFloatx80Exp(a);
2303 aSign = extractFloatx80Sign(a);
2304
2305 if (aExp == 0x7FFF) {
2306 if ((uint64_t) (aSig << 1)) {
2307 return propagateFloatx80NaNOneArg(a, status);
2308 }
2309 if (aSign) {
2310 return packFloatx80(aSign, one_exp, one_sig);
2311 }
2312 return floatx80_default_inf(0, status);
2313 }
2314
2315 if (aExp == 0 && aSig == 0) {
2316 return packFloatx80(aSign, 0, 0);
2317 }
2318
2319 user_rnd_mode = get_float_rounding_mode(status);
2320 user_rnd_prec = get_floatx80_rounding_precision(status);
2321 set_float_rounding_mode(float_round_nearest_even, status);
2322 set_floatx80_rounding_precision(floatx80_precision_x, status);
2323
2324 if (aExp >= 0x3FFD) { /* |X| >= 1/4 */
2325 compact = floatx80_make_compact(aExp, aSig);
2326
2327 if (compact <= 0x4004C215) { /* |X| <= 70 log2 */
2328 fp0 = a;
2329 fp1 = a;
2330 fp0 = floatx80_mul(fp0, float32_to_floatx80(
2331 make_float32(0x42B8AA3B), status),
2332 status); /* 64/log2 * X */
2333 n = floatx80_to_int32(fp0, status); /* int(64/log2*X) */
2334 fp0 = int32_to_floatx80(n, status);
2335
2336 j = n & 0x3F; /* J = N mod 64 */
2337 m = n / 64; /* NOTE: this is really arithmetic right shift by 6 */
2338 if (n < 0 && j) {
2339 /*
2340 * arithmetic right shift is division and
2341 * round towards minus infinity
2342 */
2343 m--;
2344 }
2345 m1 = -m;
2346 /*m += 0x3FFF; // biased exponent of 2^(M) */
2347 /*m1 += 0x3FFF; // biased exponent of -2^(-M) */
2348
2349 fp2 = fp0; /* N */
2350 fp0 = floatx80_mul(fp0, float32_to_floatx80(
2351 make_float32(0xBC317218), status),
2352 status); /* N * L1, L1 = lead(-log2/64) */
2353 l2 = packFloatx80(0, 0x3FDC, UINT64_C(0x82E308654361C4C6));
2354 fp2 = floatx80_mul(fp2, l2, status); /* N * L2, L1+L2 = -log2/64 */
2355 fp0 = floatx80_add(fp0, fp1, status); /* X + N*L1 */
2356 fp0 = floatx80_add(fp0, fp2, status); /* R */
2357
2358 fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */
2359 fp2 = float32_to_floatx80(make_float32(0x3950097B),
2360 status); /* A6 */
2361 fp2 = floatx80_mul(fp2, fp1, status); /* fp2 is S*A6 */
2362 fp3 = floatx80_mul(float32_to_floatx80(make_float32(0x3AB60B6A),
2363 status), fp1, status); /* fp3 is S*A5 */
2364 fp2 = floatx80_add(fp2, float64_to_floatx80(
2365 make_float64(0x3F81111111174385), status),
2366 status); /* fp2 IS A4+S*A6 */
2367 fp3 = floatx80_add(fp3, float64_to_floatx80(
2368 make_float64(0x3FA5555555554F5A), status),
2369 status); /* fp3 is A3+S*A5 */
2370 fp2 = floatx80_mul(fp2, fp1, status); /* fp2 IS S*(A4+S*A6) */
2371 fp3 = floatx80_mul(fp3, fp1, status); /* fp3 IS S*(A3+S*A5) */
2372 fp2 = floatx80_add(fp2, float64_to_floatx80(
2373 make_float64(0x3FC5555555555555), status),
2374 status); /* fp2 IS A2+S*(A4+S*A6) */
2375 fp3 = floatx80_add(fp3, float32_to_floatx80(
2376 make_float32(0x3F000000), status),
2377 status); /* fp3 IS A1+S*(A3+S*A5) */
2378 fp2 = floatx80_mul(fp2, fp1,
2379 status); /* fp2 IS S*(A2+S*(A4+S*A6)) */
2380 fp1 = floatx80_mul(fp1, fp3,
2381 status); /* fp1 IS S*(A1+S*(A3+S*A5)) */
2382 fp2 = floatx80_mul(fp2, fp0,
2383 status); /* fp2 IS R*S*(A2+S*(A4+S*A6)) */
2384 fp0 = floatx80_add(fp0, fp1,
2385 status); /* fp0 IS R+S*(A1+S*(A3+S*A5)) */
2386 fp0 = floatx80_add(fp0, fp2, status); /* fp0 IS EXP(R) - 1 */
2387
2388 fp0 = floatx80_mul(fp0, exp_tbl[j],
2389 status); /* 2^(J/64)*(Exp(R)-1) */
2390
2391 if (m >= 64) {
2392 fp1 = float32_to_floatx80(exp_tbl2[j], status);
2393 onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); /* -2^(-M) */
2394 fp1 = floatx80_add(fp1, onebysc, status);
2395 fp0 = floatx80_add(fp0, fp1, status);
2396 fp0 = floatx80_add(fp0, exp_tbl[j], status);
2397 } else if (m < -3) {
2398 fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j],
2399 status), status);
2400 fp0 = floatx80_add(fp0, exp_tbl[j], status);
2401 onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); /* -2^(-M) */
2402 fp0 = floatx80_add(fp0, onebysc, status);
2403 } else { /* -3 <= m <= 63 */
2404 fp1 = exp_tbl[j];
2405 fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j],
2406 status), status);
2407 onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); /* -2^(-M) */
2408 fp1 = floatx80_add(fp1, onebysc, status);
2409 fp0 = floatx80_add(fp0, fp1, status);
2410 }
2411
2412 sc = packFloatx80(0, m + 0x3FFF, one_sig);
2413
2414 set_float_rounding_mode(user_rnd_mode, status);
2415 set_floatx80_rounding_precision(user_rnd_prec, status);
2416
2417 a = floatx80_mul(fp0, sc, status);
2418
2419 float_raise(float_flag_inexact, status);
2420
2421 return a;
2422 } else { /* |X| > 70 log2 */
2423 if (aSign) {
2424 fp0 = float32_to_floatx80(make_float32(0xBF800000),
2425 status); /* -1 */
2426
2427 set_float_rounding_mode(user_rnd_mode, status);
2428 set_floatx80_rounding_precision(user_rnd_prec, status);
2429
2430 a = floatx80_add(fp0, float32_to_floatx80(
2431 make_float32(0x00800000), status),
2432 status); /* -1 + 2^(-126) */
2433
2434 float_raise(float_flag_inexact, status);
2435
2436 return a;
2437 } else {
2438 set_float_rounding_mode(user_rnd_mode, status);
2439 set_floatx80_rounding_precision(user_rnd_prec, status);
2440
2441 return floatx80_etox(a, status);
2442 }
2443 }
2444 } else { /* |X| < 1/4 */
2445 if (aExp >= 0x3FBE) {
2446 fp0 = a;
2447 fp0 = floatx80_mul(fp0, fp0, status); /* S = X*X */
2448 fp1 = float32_to_floatx80(make_float32(0x2F30CAA8),
2449 status); /* B12 */
2450 fp1 = floatx80_mul(fp1, fp0, status); /* S * B12 */
2451 fp2 = float32_to_floatx80(make_float32(0x310F8290),
2452 status); /* B11 */
2453 fp1 = floatx80_add(fp1, float32_to_floatx80(
2454 make_float32(0x32D73220), status),
2455 status); /* B10 */
2456 fp2 = floatx80_mul(fp2, fp0, status);
2457 fp1 = floatx80_mul(fp1, fp0, status);
2458 fp2 = floatx80_add(fp2, float32_to_floatx80(
2459 make_float32(0x3493F281), status),
2460 status); /* B9 */
2461 fp1 = floatx80_add(fp1, float64_to_floatx80(
2462 make_float64(0x3EC71DE3A5774682), status),
2463 status); /* B8 */
2464 fp2 = floatx80_mul(fp2, fp0, status);
2465 fp1 = floatx80_mul(fp1, fp0, status);
2466 fp2 = floatx80_add(fp2, float64_to_floatx80(
2467 make_float64(0x3EFA01A019D7CB68), status),
2468 status); /* B7 */
2469 fp1 = floatx80_add(fp1, float64_to_floatx80(
2470 make_float64(0x3F2A01A01A019DF3), status),
2471 status); /* B6 */
2472 fp2 = floatx80_mul(fp2, fp0, status);
2473 fp1 = floatx80_mul(fp1, fp0, status);
2474 fp2 = floatx80_add(fp2, float64_to_floatx80(
2475 make_float64(0x3F56C16C16C170E2), status),
2476 status); /* B5 */
2477 fp1 = floatx80_add(fp1, float64_to_floatx80(
2478 make_float64(0x3F81111111111111), status),
2479 status); /* B4 */
2480 fp2 = floatx80_mul(fp2, fp0, status);
2481 fp1 = floatx80_mul(fp1, fp0, status);
2482 fp2 = floatx80_add(fp2, float64_to_floatx80(
2483 make_float64(0x3FA5555555555555), status),
2484 status); /* B3 */
2485 fp3 = packFloatx80(0, 0x3FFC, UINT64_C(0xAAAAAAAAAAAAAAAB));
2486 fp1 = floatx80_add(fp1, fp3, status); /* B2 */
2487 fp2 = floatx80_mul(fp2, fp0, status);
2488 fp1 = floatx80_mul(fp1, fp0, status);
2489
2490 fp2 = floatx80_mul(fp2, fp0, status);
2491 fp1 = floatx80_mul(fp1, a, status);
2492
2493 fp0 = floatx80_mul(fp0, float32_to_floatx80(
2494 make_float32(0x3F000000), status),
2495 status); /* S*B1 */
2496 fp1 = floatx80_add(fp1, fp2, status); /* Q */
2497 fp0 = floatx80_add(fp0, fp1, status); /* S*B1+Q */
2498
2499 set_float_rounding_mode(user_rnd_mode, status);
2500 set_floatx80_rounding_precision(user_rnd_prec, status);
2501
2502 a = floatx80_add(fp0, a, status);
2503
2504 float_raise(float_flag_inexact, status);
2505
2506 return a;
2507 } else { /* |X| < 2^(-65) */
2508 sc = packFloatx80(1, 1, one_sig);
2509 fp0 = a;
2510
2511 if (aExp < 0x0033) { /* |X| < 2^(-16382) */
2512 fp0 = floatx80_mul(fp0, float64_to_floatx80(
2513 make_float64(0x48B0000000000000), status),
2514 status);
2515 fp0 = floatx80_add(fp0, sc, status);
2516
2517 set_float_rounding_mode(user_rnd_mode, status);
2518 set_floatx80_rounding_precision(user_rnd_prec, status);
2519
2520 a = floatx80_mul(fp0, float64_to_floatx80(
2521 make_float64(0x3730000000000000), status),
2522 status);
2523 } else {
2524 set_float_rounding_mode(user_rnd_mode, status);
2525 set_floatx80_rounding_precision(user_rnd_prec, status);
2526
2527 a = floatx80_add(fp0, sc, status);
2528 }
2529
2530 float_raise(float_flag_inexact, status);
2531
2532 return a;
2533 }
2534 }
2535 }
2536
2537 /*
2538 * Hyperbolic tangent
2539 */
2540
2541 floatx80 floatx80_tanh(floatx80 a, float_status *status)
2542 {
2543 bool aSign, vSign;
2544 int32_t aExp, vExp;
2545 uint64_t aSig, vSig;
2546
2547 FloatRoundMode user_rnd_mode;
2548 FloatX80RoundPrec user_rnd_prec;
2549
2550 int32_t compact;
2551 floatx80 fp0, fp1;
2552 uint32_t sign;
2553
2554 aSig = extractFloatx80Frac(a);
2555 aExp = extractFloatx80Exp(a);
2556 aSign = extractFloatx80Sign(a);
2557
2558 if (aExp == 0x7FFF) {
2559 if ((uint64_t) (aSig << 1)) {
2560 return propagateFloatx80NaNOneArg(a, status);
2561 }
2562 return packFloatx80(aSign, one_exp, one_sig);
2563 }
2564
2565 if (aExp == 0 && aSig == 0) {
2566 return packFloatx80(aSign, 0, 0);
2567 }
2568
2569 user_rnd_mode = get_float_rounding_mode(status);
2570 user_rnd_prec = get_floatx80_rounding_precision(status);
2571 set_float_rounding_mode(float_round_nearest_even, status);
2572 set_floatx80_rounding_precision(floatx80_precision_x, status);
2573
2574 compact = floatx80_make_compact(aExp, aSig);
2575
2576 if (compact < 0x3FD78000 || compact > 0x3FFFDDCE) {
2577 /* TANHBORS */
2578 if (compact < 0x3FFF8000) {
2579 /* TANHSM */
2580 set_float_rounding_mode(user_rnd_mode, status);
2581 set_floatx80_rounding_precision(user_rnd_prec, status);
2582
2583 a = floatx80_move(a, status);
2584
2585 float_raise(float_flag_inexact, status);
2586
2587 return a;
2588 } else {
2589 if (compact > 0x40048AA1) {
2590 /* TANHHUGE */
2591 sign = 0x3F800000;
2592 sign |= aSign ? 0x80000000 : 0x00000000;
2593 fp0 = float32_to_floatx80(make_float32(sign), status);
2594 sign &= 0x80000000;
2595 sign ^= 0x80800000; /* -SIGN(X)*EPS */
2596
2597 set_float_rounding_mode(user_rnd_mode, status);
2598 set_floatx80_rounding_precision(user_rnd_prec, status);
2599
2600 a = floatx80_add(fp0, float32_to_floatx80(make_float32(sign),
2601 status), status);
2602
2603 float_raise(float_flag_inexact, status);
2604
2605 return a;
2606 } else {
2607 fp0 = packFloatx80(0, aExp + 1, aSig); /* Y = 2|X| */
2608 fp0 = floatx80_etox(fp0, status); /* FP0 IS EXP(Y) */
2609 fp0 = floatx80_add(fp0, float32_to_floatx80(
2610 make_float32(0x3F800000),
2611 status), status); /* EXP(Y)+1 */
2612 sign = aSign ? 0x80000000 : 0x00000000;
2613 fp1 = floatx80_div(float32_to_floatx80(make_float32(
2614 sign ^ 0xC0000000), status), fp0,
2615 status); /* -SIGN(X)*2 / [EXP(Y)+1] */
2616 fp0 = float32_to_floatx80(make_float32(sign | 0x3F800000),
2617 status); /* SIGN */
2618
2619 set_float_rounding_mode(user_rnd_mode, status);
2620 set_floatx80_rounding_precision(user_rnd_prec, status);
2621
2622 a = floatx80_add(fp1, fp0, status);
2623
2624 float_raise(float_flag_inexact, status);
2625
2626 return a;
2627 }
2628 }
2629 } else { /* 2**(-40) < |X| < (5/2)LOG2 */
2630 fp0 = packFloatx80(0, aExp + 1, aSig); /* Y = 2|X| */
2631 fp0 = floatx80_etoxm1(fp0, status); /* FP0 IS Z = EXPM1(Y) */
2632 fp1 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x40000000),
2633 status),
2634 status); /* Z+2 */
2635
2636 vSign = extractFloatx80Sign(fp1);
2637 vExp = extractFloatx80Exp(fp1);
2638 vSig = extractFloatx80Frac(fp1);
2639
2640 fp1 = packFloatx80(vSign ^ aSign, vExp, vSig);
2641
2642 set_float_rounding_mode(user_rnd_mode, status);
2643 set_floatx80_rounding_precision(user_rnd_prec, status);
2644
2645 a = floatx80_div(fp0, fp1, status);
2646
2647 float_raise(float_flag_inexact, status);
2648
2649 return a;
2650 }
2651 }
2652
2653 /*
2654 * Hyperbolic sine
2655 */
2656
2657 floatx80 floatx80_sinh(floatx80 a, float_status *status)
2658 {
2659 bool aSign;
2660 int32_t aExp;
2661 uint64_t aSig;
2662
2663 FloatRoundMode user_rnd_mode;
2664 FloatX80RoundPrec user_rnd_prec;
2665
2666 int32_t compact;
2667 floatx80 fp0, fp1, fp2;
2668 float32 fact;
2669
2670 aSig = extractFloatx80Frac(a);
2671 aExp = extractFloatx80Exp(a);
2672 aSign = extractFloatx80Sign(a);
2673
2674 if (aExp == 0x7FFF) {
2675 if ((uint64_t) (aSig << 1)) {
2676 return propagateFloatx80NaNOneArg(a, status);
2677 }
2678 return floatx80_default_inf(aSign, status);
2679 }
2680
2681 if (aExp == 0 && aSig == 0) {
2682 return packFloatx80(aSign, 0, 0);
2683 }
2684
2685 user_rnd_mode = get_float_rounding_mode(status);
2686 user_rnd_prec = get_floatx80_rounding_precision(status);
2687 set_float_rounding_mode(float_round_nearest_even, status);
2688 set_floatx80_rounding_precision(floatx80_precision_x, status);
2689
2690 compact = floatx80_make_compact(aExp, aSig);
2691
2692 if (compact > 0x400CB167) {
2693 /* SINHBIG */
2694 if (compact > 0x400CB2B3) {
2695 set_float_rounding_mode(user_rnd_mode, status);
2696 set_floatx80_rounding_precision(user_rnd_prec, status);
2697
2698 return roundAndPackFloatx80(get_floatx80_rounding_precision(status),
2699 aSign, 0x8000, aSig, 0, status);
2700 } else {
2701 fp0 = floatx80_abs(a); /* Y = |X| */
2702 fp0 = floatx80_sub(fp0, float64_to_floatx80(
2703 make_float64(0x40C62D38D3D64634), status),
2704 status); /* (|X|-16381LOG2_LEAD) */
2705 fp0 = floatx80_sub(fp0, float64_to_floatx80(
2706 make_float64(0x3D6F90AEB1E75CC7), status),
2707 status); /* |X| - 16381 LOG2, ACCURATE */
2708 fp0 = floatx80_etox(fp0, status);
2709 fp2 = packFloatx80(aSign, 0x7FFB, one_sig);
2710
2711 set_float_rounding_mode(user_rnd_mode, status);
2712 set_floatx80_rounding_precision(user_rnd_prec, status);
2713
2714 a = floatx80_mul(fp0, fp2, status);
2715
2716 float_raise(float_flag_inexact, status);
2717
2718 return a;
2719 }
2720 } else { /* |X| < 16380 LOG2 */
2721 fp0 = floatx80_abs(a); /* Y = |X| */
2722 fp0 = floatx80_etoxm1(fp0, status); /* FP0 IS Z = EXPM1(Y) */
2723 fp1 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000),
2724 status), status); /* 1+Z */
2725 fp2 = fp0;
2726 fp0 = floatx80_div(fp0, fp1, status); /* Z/(1+Z) */
2727 fp0 = floatx80_add(fp0, fp2, status);
2728
2729 fact = packFloat32(aSign, 0x7E, 0);
2730
2731 set_float_rounding_mode(user_rnd_mode, status);
2732 set_floatx80_rounding_precision(user_rnd_prec, status);
2733
2734 a = floatx80_mul(fp0, float32_to_floatx80(fact, status), status);
2735
2736 float_raise(float_flag_inexact, status);
2737
2738 return a;
2739 }
2740 }
2741
2742 /*
2743 * Hyperbolic cosine
2744 */
2745
2746 floatx80 floatx80_cosh(floatx80 a, float_status *status)
2747 {
2748 int32_t aExp;
2749 uint64_t aSig;
2750
2751 FloatRoundMode user_rnd_mode;
2752 FloatX80RoundPrec user_rnd_prec;
2753
2754 int32_t compact;
2755 floatx80 fp0, fp1;
2756
2757 aSig = extractFloatx80Frac(a);
2758 aExp = extractFloatx80Exp(a);
2759
2760 if (aExp == 0x7FFF) {
2761 if ((uint64_t) (aSig << 1)) {
2762 return propagateFloatx80NaNOneArg(a, status);
2763 }
2764 return floatx80_default_inf(0, status);
2765 }
2766
2767 if (aExp == 0 && aSig == 0) {
2768 return packFloatx80(0, one_exp, one_sig);
2769 }
2770
2771 user_rnd_mode = get_float_rounding_mode(status);
2772 user_rnd_prec = get_floatx80_rounding_precision(status);
2773 set_float_rounding_mode(float_round_nearest_even, status);
2774 set_floatx80_rounding_precision(floatx80_precision_x, status);
2775
2776 compact = floatx80_make_compact(aExp, aSig);
2777
2778 if (compact > 0x400CB167) {
2779 if (compact > 0x400CB2B3) {
2780 set_float_rounding_mode(user_rnd_mode, status);
2781 set_floatx80_rounding_precision(user_rnd_prec, status);
2782 return roundAndPackFloatx80(get_floatx80_rounding_precision(status),
2783 0,
2784 0x8000, one_sig, 0, status);
2785 } else {
2786 fp0 = packFloatx80(0, aExp, aSig);
2787 fp0 = floatx80_sub(fp0, float64_to_floatx80(
2788 make_float64(0x40C62D38D3D64634), status),
2789 status);
2790 fp0 = floatx80_sub(fp0, float64_to_floatx80(
2791 make_float64(0x3D6F90AEB1E75CC7), status),
2792 status);
2793 fp0 = floatx80_etox(fp0, status);
2794 fp1 = packFloatx80(0, 0x7FFB, one_sig);
2795
2796 set_float_rounding_mode(user_rnd_mode, status);
2797 set_floatx80_rounding_precision(user_rnd_prec, status);
2798
2799 a = floatx80_mul(fp0, fp1, status);
2800
2801 float_raise(float_flag_inexact, status);
2802
2803 return a;
2804 }
2805 }
2806
2807 fp0 = packFloatx80(0, aExp, aSig); /* |X| */
2808 fp0 = floatx80_etox(fp0, status); /* EXP(|X|) */
2809 fp0 = floatx80_mul(fp0, float32_to_floatx80(make_float32(0x3F000000),
2810 status), status); /* (1/2)*EXP(|X|) */
2811 fp1 = float32_to_floatx80(make_float32(0x3E800000), status); /* 1/4 */
2812 fp1 = floatx80_div(fp1, fp0, status); /* 1/(2*EXP(|X|)) */
2813
2814 set_float_rounding_mode(user_rnd_mode, status);
2815 set_floatx80_rounding_precision(user_rnd_prec, status);
2816
2817 a = floatx80_add(fp0, fp1, status);
2818
2819 float_raise(float_flag_inexact, status);
2820
2821 return a;
2822 }