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
propagateFloatx80NaNOneArg(floatx80 a,float_status * status)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 (status->default_nan_mode) {
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
floatx80_getman(floatx80 a,float_status * status)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(status->floatx80_rounding_precision, aSign,
76 0x3FFF, aSig, 0, status);
77 }
78
79 /*
80 * Returns the exponent of the extended double-precision floating-point
81 * value `a' as an extended double-precision value.
82 */
83
floatx80_getexp(floatx80 a,float_status * status)84 floatx80 floatx80_getexp(floatx80 a, float_status *status)
85 {
86 bool aSign;
87 int32_t aExp;
88 uint64_t aSig;
89
90 aSig = extractFloatx80Frac(a);
91 aExp = extractFloatx80Exp(a);
92 aSign = extractFloatx80Sign(a);
93
94 if (aExp == 0x7FFF) {
95 if ((uint64_t) (aSig << 1)) {
96 return propagateFloatx80NaNOneArg(a , status);
97 }
98 float_raise(float_flag_invalid , status);
99 return floatx80_default_nan(status);
100 }
101
102 if (aExp == 0) {
103 if (aSig == 0) {
104 return packFloatx80(aSign, 0, 0);
105 }
106 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
107 }
108
109 return int32_to_floatx80(aExp - 0x3FFF, status);
110 }
111
112 /*
113 * Scales extended double-precision floating-point value in operand `a' by
114 * value `b'. The function truncates the value in the second operand 'b' to
115 * an integral value and adds that value to the exponent of the operand 'a'.
116 * The operation performed according to the IEC/IEEE Standard for Binary
117 * Floating-Point Arithmetic.
118 */
119
floatx80_scale(floatx80 a,floatx80 b,float_status * status)120 floatx80 floatx80_scale(floatx80 a, floatx80 b, float_status *status)
121 {
122 bool aSign, bSign;
123 int32_t aExp, bExp, shiftCount;
124 uint64_t aSig, bSig;
125
126 aSig = extractFloatx80Frac(a);
127 aExp = extractFloatx80Exp(a);
128 aSign = extractFloatx80Sign(a);
129 bSig = extractFloatx80Frac(b);
130 bExp = extractFloatx80Exp(b);
131 bSign = extractFloatx80Sign(b);
132
133 if (bExp == 0x7FFF) {
134 if ((uint64_t) (bSig << 1) ||
135 ((aExp == 0x7FFF) && (uint64_t) (aSig << 1))) {
136 return propagateFloatx80NaN(a, b, status);
137 }
138 float_raise(float_flag_invalid , status);
139 return floatx80_default_nan(status);
140 }
141 if (aExp == 0x7FFF) {
142 if ((uint64_t) (aSig << 1)) {
143 return propagateFloatx80NaN(a, b, status);
144 }
145 return floatx80_default_inf(aSign, status);
146 }
147 if (aExp == 0) {
148 if (aSig == 0) {
149 return packFloatx80(aSign, 0, 0);
150 }
151 if (bExp < 0x3FFF) {
152 return a;
153 }
154 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
155 }
156
157 if (bExp < 0x3FFF) {
158 return a;
159 }
160
161 if (0x400F < bExp) {
162 aExp = bSign ? -0x6001 : 0xE000;
163 return roundAndPackFloatx80(status->floatx80_rounding_precision,
164 aSign, aExp, aSig, 0, status);
165 }
166
167 shiftCount = 0x403E - bExp;
168 bSig >>= shiftCount;
169 aExp = bSign ? (aExp - bSig) : (aExp + bSig);
170
171 return roundAndPackFloatx80(status->floatx80_rounding_precision,
172 aSign, aExp, aSig, 0, status);
173 }
174
floatx80_move(floatx80 a,float_status * status)175 floatx80 floatx80_move(floatx80 a, float_status *status)
176 {
177 bool aSign;
178 int32_t aExp;
179 uint64_t aSig;
180
181 aSig = extractFloatx80Frac(a);
182 aExp = extractFloatx80Exp(a);
183 aSign = extractFloatx80Sign(a);
184
185 if (aExp == 0x7FFF) {
186 if ((uint64_t)(aSig << 1)) {
187 return propagateFloatx80NaNOneArg(a, status);
188 }
189 return a;
190 }
191 if (aExp == 0) {
192 if (aSig == 0) {
193 return a;
194 }
195 normalizeRoundAndPackFloatx80(status->floatx80_rounding_precision,
196 aSign, aExp, aSig, 0, status);
197 }
198 return roundAndPackFloatx80(status->floatx80_rounding_precision, aSign,
199 aExp, aSig, 0, status);
200 }
201
202 /*
203 * Algorithms for transcendental functions supported by MC68881 and MC68882
204 * mathematical coprocessors. The functions are derived from FPSP library.
205 */
206
207 #define one_exp 0x3FFF
208 #define one_sig UINT64_C(0x8000000000000000)
209
210 /*
211 * Function for compactifying extended double-precision floating point values.
212 */
213
floatx80_make_compact(int32_t aExp,uint64_t aSig)214 static int32_t floatx80_make_compact(int32_t aExp, uint64_t aSig)
215 {
216 return (aExp << 16) | (aSig >> 48);
217 }
218
219 /*
220 * Log base e of x plus 1
221 */
222
floatx80_lognp1(floatx80 a,float_status * status)223 floatx80 floatx80_lognp1(floatx80 a, float_status *status)
224 {
225 bool aSign;
226 int32_t aExp;
227 uint64_t aSig, fSig;
228
229 FloatRoundMode user_rnd_mode;
230 FloatX80RoundPrec user_rnd_prec;
231
232 int32_t compact, j, k;
233 floatx80 fp0, fp1, fp2, fp3, f, logof2, klog2, saveu;
234
235 aSig = extractFloatx80Frac(a);
236 aExp = extractFloatx80Exp(a);
237 aSign = extractFloatx80Sign(a);
238
239 if (aExp == 0x7FFF) {
240 if ((uint64_t) (aSig << 1)) {
241 propagateFloatx80NaNOneArg(a, status);
242 }
243 if (aSign) {
244 float_raise(float_flag_invalid, status);
245 return floatx80_default_nan(status);
246 }
247 return floatx80_default_inf(0, status);
248 }
249
250 if (aExp == 0 && aSig == 0) {
251 return packFloatx80(aSign, 0, 0);
252 }
253
254 if (aSign && aExp >= one_exp) {
255 if (aExp == one_exp && aSig == one_sig) {
256 float_raise(float_flag_divbyzero, status);
257 return floatx80_default_inf(aSign, status);
258 }
259 float_raise(float_flag_invalid, status);
260 return floatx80_default_nan(status);
261 }
262
263 if (aExp < 0x3f99 || (aExp == 0x3f99 && aSig == one_sig)) {
264 /* <= min threshold */
265 float_raise(float_flag_inexact, status);
266 return floatx80_move(a, status);
267 }
268
269 user_rnd_mode = status->float_rounding_mode;
270 user_rnd_prec = status->floatx80_rounding_precision;
271 status->float_rounding_mode = float_round_nearest_even;
272 status->floatx80_rounding_precision = floatx80_precision_x;
273
274 compact = floatx80_make_compact(aExp, aSig);
275
276 fp0 = a; /* Z */
277 fp1 = a;
278
279 fp0 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000),
280 status), status); /* X = (1+Z) */
281
282 aExp = extractFloatx80Exp(fp0);
283 aSig = extractFloatx80Frac(fp0);
284
285 compact = floatx80_make_compact(aExp, aSig);
286
287 if (compact < 0x3FFE8000 || compact > 0x3FFFC000) {
288 /* |X| < 1/2 or |X| > 3/2 */
289 k = aExp - 0x3FFF;
290 fp1 = int32_to_floatx80(k, status);
291
292 fSig = (aSig & UINT64_C(0xFE00000000000000)) | UINT64_C(0x0100000000000000);
293 j = (fSig >> 56) & 0x7E; /* DISPLACEMENT FOR 1/F */
294
295 f = packFloatx80(0, 0x3FFF, fSig); /* F */
296 fp0 = packFloatx80(0, 0x3FFF, aSig); /* Y */
297
298 fp0 = floatx80_sub(fp0, f, status); /* Y-F */
299
300 lp1cont1:
301 /* LP1CONT1 */
302 fp0 = floatx80_mul(fp0, log_tbl[j], status); /* FP0 IS U = (Y-F)/F */
303 logof2 = packFloatx80(0, 0x3FFE, UINT64_C(0xB17217F7D1CF79AC));
304 klog2 = floatx80_mul(fp1, logof2, status); /* FP1 IS K*LOG2 */
305 fp2 = floatx80_mul(fp0, fp0, status); /* FP2 IS V=U*U */
306
307 fp3 = fp2;
308 fp1 = fp2;
309
310 fp1 = floatx80_mul(fp1, float64_to_floatx80(
311 make_float64(0x3FC2499AB5E4040B), status),
312 status); /* V*A6 */
313 fp2 = floatx80_mul(fp2, float64_to_floatx80(
314 make_float64(0xBFC555B5848CB7DB), status),
315 status); /* V*A5 */
316 fp1 = floatx80_add(fp1, float64_to_floatx80(
317 make_float64(0x3FC99999987D8730), status),
318 status); /* A4+V*A6 */
319 fp2 = floatx80_add(fp2, float64_to_floatx80(
320 make_float64(0xBFCFFFFFFF6F7E97), status),
321 status); /* A3+V*A5 */
322 fp1 = floatx80_mul(fp1, fp3, status); /* V*(A4+V*A6) */
323 fp2 = floatx80_mul(fp2, fp3, status); /* V*(A3+V*A5) */
324 fp1 = floatx80_add(fp1, float64_to_floatx80(
325 make_float64(0x3FD55555555555A4), status),
326 status); /* A2+V*(A4+V*A6) */
327 fp2 = floatx80_add(fp2, float64_to_floatx80(
328 make_float64(0xBFE0000000000008), status),
329 status); /* A1+V*(A3+V*A5) */
330 fp1 = floatx80_mul(fp1, fp3, status); /* V*(A2+V*(A4+V*A6)) */
331 fp2 = floatx80_mul(fp2, fp3, status); /* V*(A1+V*(A3+V*A5)) */
332 fp1 = floatx80_mul(fp1, fp0, status); /* U*V*(A2+V*(A4+V*A6)) */
333 fp0 = floatx80_add(fp0, fp2, status); /* U+V*(A1+V*(A3+V*A5)) */
334
335 fp1 = floatx80_add(fp1, log_tbl[j + 1],
336 status); /* LOG(F)+U*V*(A2+V*(A4+V*A6)) */
337 fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS LOG(F) + LOG(1+U) */
338
339 status->float_rounding_mode = user_rnd_mode;
340 status->floatx80_rounding_precision = user_rnd_prec;
341
342 a = floatx80_add(fp0, klog2, status);
343
344 float_raise(float_flag_inexact, status);
345
346 return a;
347 } else if (compact < 0x3FFEF07D || compact > 0x3FFF8841) {
348 /* |X| < 1/16 or |X| > -1/16 */
349 /* LP1CARE */
350 fSig = (aSig & UINT64_C(0xFE00000000000000)) | UINT64_C(0x0100000000000000);
351 f = packFloatx80(0, 0x3FFF, fSig); /* F */
352 j = (fSig >> 56) & 0x7E; /* DISPLACEMENT FOR 1/F */
353
354 if (compact >= 0x3FFF8000) { /* 1+Z >= 1 */
355 /* KISZERO */
356 fp0 = floatx80_sub(float32_to_floatx80(make_float32(0x3F800000),
357 status), f, status); /* 1-F */
358 fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS Y-F = (1-F)+Z */
359 fp1 = packFloatx80(0, 0, 0); /* K = 0 */
360 } else {
361 /* KISNEG */
362 fp0 = floatx80_sub(float32_to_floatx80(make_float32(0x40000000),
363 status), f, status); /* 2-F */
364 fp1 = floatx80_add(fp1, fp1, status); /* 2Z */
365 fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS Y-F = (2-F)+2Z */
366 fp1 = packFloatx80(1, one_exp, one_sig); /* K = -1 */
367 }
368 goto lp1cont1;
369 } else {
370 /* LP1ONE16 */
371 fp1 = floatx80_add(fp1, fp1, status); /* FP1 IS 2Z */
372 fp0 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000),
373 status), status); /* FP0 IS 1+X */
374
375 /* LP1CONT2 */
376 fp1 = floatx80_div(fp1, fp0, status); /* U */
377 saveu = fp1;
378 fp0 = floatx80_mul(fp1, fp1, status); /* FP0 IS V = U*U */
379 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS W = V*V */
380
381 fp3 = float64_to_floatx80(make_float64(0x3F175496ADD7DAD6),
382 status); /* B5 */
383 fp2 = float64_to_floatx80(make_float64(0x3F3C71C2FE80C7E0),
384 status); /* B4 */
385 fp3 = floatx80_mul(fp3, fp1, status); /* W*B5 */
386 fp2 = floatx80_mul(fp2, fp1, status); /* W*B4 */
387 fp3 = floatx80_add(fp3, float64_to_floatx80(
388 make_float64(0x3F624924928BCCFF), status),
389 status); /* B3+W*B5 */
390 fp2 = floatx80_add(fp2, float64_to_floatx80(
391 make_float64(0x3F899999999995EC), status),
392 status); /* B2+W*B4 */
393 fp1 = floatx80_mul(fp1, fp3, status); /* W*(B3+W*B5) */
394 fp2 = floatx80_mul(fp2, fp0, status); /* V*(B2+W*B4) */
395 fp1 = floatx80_add(fp1, float64_to_floatx80(
396 make_float64(0x3FB5555555555555), status),
397 status); /* B1+W*(B3+W*B5) */
398
399 fp0 = floatx80_mul(fp0, saveu, status); /* FP0 IS U*V */
400 fp1 = floatx80_add(fp1, fp2,
401 status); /* B1+W*(B3+W*B5) + V*(B2+W*B4) */
402 fp0 = floatx80_mul(fp0, fp1,
403 status); /* U*V*([B1+W*(B3+W*B5)] + [V*(B2+W*B4)]) */
404
405 status->float_rounding_mode = user_rnd_mode;
406 status->floatx80_rounding_precision = user_rnd_prec;
407
408 a = floatx80_add(fp0, saveu, status);
409
410 /*if (!floatx80_is_zero(a)) { */
411 float_raise(float_flag_inexact, status);
412 /*} */
413
414 return a;
415 }
416 }
417
418 /*
419 * Log base e
420 */
421
floatx80_logn(floatx80 a,float_status * status)422 floatx80 floatx80_logn(floatx80 a, float_status *status)
423 {
424 bool aSign;
425 int32_t aExp;
426 uint64_t aSig, fSig;
427
428 FloatRoundMode user_rnd_mode;
429 FloatX80RoundPrec user_rnd_prec;
430
431 int32_t compact, j, k, adjk;
432 floatx80 fp0, fp1, fp2, fp3, f, logof2, klog2, saveu;
433
434 aSig = extractFloatx80Frac(a);
435 aExp = extractFloatx80Exp(a);
436 aSign = extractFloatx80Sign(a);
437
438 if (aExp == 0x7FFF) {
439 if ((uint64_t) (aSig << 1)) {
440 propagateFloatx80NaNOneArg(a, status);
441 }
442 if (aSign == 0) {
443 return floatx80_default_inf(0, status);
444 }
445 }
446
447 adjk = 0;
448
449 if (aExp == 0) {
450 if (aSig == 0) { /* zero */
451 float_raise(float_flag_divbyzero, status);
452 return floatx80_default_inf(1, status);
453 }
454 if ((aSig & one_sig) == 0) { /* denormal */
455 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
456 adjk = -100;
457 aExp += 100;
458 a = packFloatx80(aSign, aExp, aSig);
459 }
460 }
461
462 if (aSign) {
463 float_raise(float_flag_invalid, status);
464 return floatx80_default_nan(status);
465 }
466
467 user_rnd_mode = status->float_rounding_mode;
468 user_rnd_prec = status->floatx80_rounding_precision;
469 status->float_rounding_mode = float_round_nearest_even;
470 status->floatx80_rounding_precision = floatx80_precision_x;
471
472 compact = floatx80_make_compact(aExp, aSig);
473
474 if (compact < 0x3FFEF07D || compact > 0x3FFF8841) {
475 /* |X| < 15/16 or |X| > 17/16 */
476 k = aExp - 0x3FFF;
477 k += adjk;
478 fp1 = int32_to_floatx80(k, status);
479
480 fSig = (aSig & UINT64_C(0xFE00000000000000)) | UINT64_C(0x0100000000000000);
481 j = (fSig >> 56) & 0x7E; /* DISPLACEMENT FOR 1/F */
482
483 f = packFloatx80(0, 0x3FFF, fSig); /* F */
484 fp0 = packFloatx80(0, 0x3FFF, aSig); /* Y */
485
486 fp0 = floatx80_sub(fp0, f, status); /* Y-F */
487
488 /* LP1CONT1 */
489 fp0 = floatx80_mul(fp0, log_tbl[j], status); /* FP0 IS U = (Y-F)/F */
490 logof2 = packFloatx80(0, 0x3FFE, UINT64_C(0xB17217F7D1CF79AC));
491 klog2 = floatx80_mul(fp1, logof2, status); /* FP1 IS K*LOG2 */
492 fp2 = floatx80_mul(fp0, fp0, status); /* FP2 IS V=U*U */
493
494 fp3 = fp2;
495 fp1 = fp2;
496
497 fp1 = floatx80_mul(fp1, float64_to_floatx80(
498 make_float64(0x3FC2499AB5E4040B), status),
499 status); /* V*A6 */
500 fp2 = floatx80_mul(fp2, float64_to_floatx80(
501 make_float64(0xBFC555B5848CB7DB), status),
502 status); /* V*A5 */
503 fp1 = floatx80_add(fp1, float64_to_floatx80(
504 make_float64(0x3FC99999987D8730), status),
505 status); /* A4+V*A6 */
506 fp2 = floatx80_add(fp2, float64_to_floatx80(
507 make_float64(0xBFCFFFFFFF6F7E97), status),
508 status); /* A3+V*A5 */
509 fp1 = floatx80_mul(fp1, fp3, status); /* V*(A4+V*A6) */
510 fp2 = floatx80_mul(fp2, fp3, status); /* V*(A3+V*A5) */
511 fp1 = floatx80_add(fp1, float64_to_floatx80(
512 make_float64(0x3FD55555555555A4), status),
513 status); /* A2+V*(A4+V*A6) */
514 fp2 = floatx80_add(fp2, float64_to_floatx80(
515 make_float64(0xBFE0000000000008), status),
516 status); /* A1+V*(A3+V*A5) */
517 fp1 = floatx80_mul(fp1, fp3, status); /* V*(A2+V*(A4+V*A6)) */
518 fp2 = floatx80_mul(fp2, fp3, status); /* V*(A1+V*(A3+V*A5)) */
519 fp1 = floatx80_mul(fp1, fp0, status); /* U*V*(A2+V*(A4+V*A6)) */
520 fp0 = floatx80_add(fp0, fp2, status); /* U+V*(A1+V*(A3+V*A5)) */
521
522 fp1 = floatx80_add(fp1, log_tbl[j + 1],
523 status); /* LOG(F)+U*V*(A2+V*(A4+V*A6)) */
524 fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS LOG(F) + LOG(1+U) */
525
526 status->float_rounding_mode = user_rnd_mode;
527 status->floatx80_rounding_precision = user_rnd_prec;
528
529 a = floatx80_add(fp0, klog2, status);
530
531 float_raise(float_flag_inexact, status);
532
533 return a;
534 } else { /* |X-1| >= 1/16 */
535 fp0 = a;
536 fp1 = a;
537 fp1 = floatx80_sub(fp1, float32_to_floatx80(make_float32(0x3F800000),
538 status), status); /* FP1 IS X-1 */
539 fp0 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000),
540 status), status); /* FP0 IS X+1 */
541 fp1 = floatx80_add(fp1, fp1, status); /* FP1 IS 2(X-1) */
542
543 /* LP1CONT2 */
544 fp1 = floatx80_div(fp1, fp0, status); /* U */
545 saveu = fp1;
546 fp0 = floatx80_mul(fp1, fp1, status); /* FP0 IS V = U*U */
547 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS W = V*V */
548
549 fp3 = float64_to_floatx80(make_float64(0x3F175496ADD7DAD6),
550 status); /* B5 */
551 fp2 = float64_to_floatx80(make_float64(0x3F3C71C2FE80C7E0),
552 status); /* B4 */
553 fp3 = floatx80_mul(fp3, fp1, status); /* W*B5 */
554 fp2 = floatx80_mul(fp2, fp1, status); /* W*B4 */
555 fp3 = floatx80_add(fp3, float64_to_floatx80(
556 make_float64(0x3F624924928BCCFF), status),
557 status); /* B3+W*B5 */
558 fp2 = floatx80_add(fp2, float64_to_floatx80(
559 make_float64(0x3F899999999995EC), status),
560 status); /* B2+W*B4 */
561 fp1 = floatx80_mul(fp1, fp3, status); /* W*(B3+W*B5) */
562 fp2 = floatx80_mul(fp2, fp0, status); /* V*(B2+W*B4) */
563 fp1 = floatx80_add(fp1, float64_to_floatx80(
564 make_float64(0x3FB5555555555555), status),
565 status); /* B1+W*(B3+W*B5) */
566
567 fp0 = floatx80_mul(fp0, saveu, status); /* FP0 IS U*V */
568 fp1 = floatx80_add(fp1, fp2, status); /* B1+W*(B3+W*B5) + V*(B2+W*B4) */
569 fp0 = floatx80_mul(fp0, fp1,
570 status); /* U*V*([B1+W*(B3+W*B5)] + [V*(B2+W*B4)]) */
571
572 status->float_rounding_mode = user_rnd_mode;
573 status->floatx80_rounding_precision = user_rnd_prec;
574
575 a = floatx80_add(fp0, saveu, status);
576
577 /*if (!floatx80_is_zero(a)) { */
578 float_raise(float_flag_inexact, status);
579 /*} */
580
581 return a;
582 }
583 }
584
585 /*
586 * Log base 10
587 */
588
floatx80_log10(floatx80 a,float_status * status)589 floatx80 floatx80_log10(floatx80 a, float_status *status)
590 {
591 bool aSign;
592 int32_t aExp;
593 uint64_t aSig;
594
595 FloatRoundMode user_rnd_mode;
596 FloatX80RoundPrec user_rnd_prec;
597
598 floatx80 fp0, fp1;
599
600 aSig = extractFloatx80Frac(a);
601 aExp = extractFloatx80Exp(a);
602 aSign = extractFloatx80Sign(a);
603
604 if (aExp == 0x7FFF) {
605 if ((uint64_t) (aSig << 1)) {
606 propagateFloatx80NaNOneArg(a, status);
607 }
608 if (aSign == 0) {
609 return floatx80_default_inf(0, status);
610 }
611 }
612
613 if (aExp == 0 && aSig == 0) {
614 float_raise(float_flag_divbyzero, status);
615 return floatx80_default_inf(1, status);
616 }
617
618 if (aSign) {
619 float_raise(float_flag_invalid, status);
620 return floatx80_default_nan(status);
621 }
622
623 user_rnd_mode = status->float_rounding_mode;
624 user_rnd_prec = status->floatx80_rounding_precision;
625 status->float_rounding_mode = float_round_nearest_even;
626 status->floatx80_rounding_precision = floatx80_precision_x;
627
628 fp0 = floatx80_logn(a, status);
629 fp1 = packFloatx80(0, 0x3FFD, UINT64_C(0xDE5BD8A937287195)); /* INV_L10 */
630
631 status->float_rounding_mode = user_rnd_mode;
632 status->floatx80_rounding_precision = user_rnd_prec;
633
634 a = floatx80_mul(fp0, fp1, status); /* LOGN(X)*INV_L10 */
635
636 float_raise(float_flag_inexact, status);
637
638 return a;
639 }
640
641 /*
642 * Log base 2
643 */
644
floatx80_log2(floatx80 a,float_status * status)645 floatx80 floatx80_log2(floatx80 a, float_status *status)
646 {
647 bool aSign;
648 int32_t aExp;
649 uint64_t aSig;
650
651 FloatRoundMode user_rnd_mode;
652 FloatX80RoundPrec user_rnd_prec;
653
654 floatx80 fp0, fp1;
655
656 aSig = extractFloatx80Frac(a);
657 aExp = extractFloatx80Exp(a);
658 aSign = extractFloatx80Sign(a);
659
660 if (aExp == 0x7FFF) {
661 if ((uint64_t) (aSig << 1)) {
662 propagateFloatx80NaNOneArg(a, status);
663 }
664 if (aSign == 0) {
665 return floatx80_default_inf(0, status);
666 }
667 }
668
669 if (aExp == 0) {
670 if (aSig == 0) {
671 float_raise(float_flag_divbyzero, status);
672 return floatx80_default_inf(1, status);
673 }
674 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
675 }
676
677 if (aSign) {
678 float_raise(float_flag_invalid, status);
679 return floatx80_default_nan(status);
680 }
681
682 user_rnd_mode = status->float_rounding_mode;
683 user_rnd_prec = status->floatx80_rounding_precision;
684 status->float_rounding_mode = float_round_nearest_even;
685 status->floatx80_rounding_precision = floatx80_precision_x;
686
687 if (aSig == one_sig) { /* X is 2^k */
688 status->float_rounding_mode = user_rnd_mode;
689 status->floatx80_rounding_precision = user_rnd_prec;
690
691 a = int32_to_floatx80(aExp - 0x3FFF, status);
692 } else {
693 fp0 = floatx80_logn(a, status);
694 fp1 = packFloatx80(0, 0x3FFF, UINT64_C(0xB8AA3B295C17F0BC)); /* INV_L2 */
695
696 status->float_rounding_mode = user_rnd_mode;
697 status->floatx80_rounding_precision = user_rnd_prec;
698
699 a = floatx80_mul(fp0, fp1, status); /* LOGN(X)*INV_L2 */
700 }
701
702 float_raise(float_flag_inexact, status);
703
704 return a;
705 }
706
707 /*
708 * e to x
709 */
710
floatx80_etox(floatx80 a,float_status * status)711 floatx80 floatx80_etox(floatx80 a, float_status *status)
712 {
713 bool aSign;
714 int32_t aExp;
715 uint64_t aSig;
716
717 FloatRoundMode user_rnd_mode;
718 FloatX80RoundPrec user_rnd_prec;
719
720 int32_t compact, n, j, k, m, m1;
721 floatx80 fp0, fp1, fp2, fp3, l2, scale, adjscale;
722 bool adjflag;
723
724 aSig = extractFloatx80Frac(a);
725 aExp = extractFloatx80Exp(a);
726 aSign = extractFloatx80Sign(a);
727
728 if (aExp == 0x7FFF) {
729 if ((uint64_t) (aSig << 1)) {
730 return propagateFloatx80NaNOneArg(a, status);
731 }
732 if (aSign) {
733 return packFloatx80(0, 0, 0);
734 }
735 return floatx80_default_inf(0, status);
736 }
737
738 if (aExp == 0 && aSig == 0) {
739 return packFloatx80(0, one_exp, one_sig);
740 }
741
742 user_rnd_mode = status->float_rounding_mode;
743 user_rnd_prec = status->floatx80_rounding_precision;
744 status->float_rounding_mode = float_round_nearest_even;
745 status->floatx80_rounding_precision = floatx80_precision_x;
746
747 adjflag = 0;
748
749 if (aExp >= 0x3FBE) { /* |X| >= 2^(-65) */
750 compact = floatx80_make_compact(aExp, aSig);
751
752 if (compact < 0x400CB167) { /* |X| < 16380 log2 */
753 fp0 = a;
754 fp1 = a;
755 fp0 = floatx80_mul(fp0, float32_to_floatx80(
756 make_float32(0x42B8AA3B), status),
757 status); /* 64/log2 * X */
758 adjflag = 0;
759 n = floatx80_to_int32(fp0, status); /* int(64/log2*X) */
760 fp0 = int32_to_floatx80(n, status);
761
762 j = n & 0x3F; /* J = N mod 64 */
763 m = n / 64; /* NOTE: this is really arithmetic right shift by 6 */
764 if (n < 0 && j) {
765 /*
766 * arithmetic right shift is division and
767 * round towards minus infinity
768 */
769 m--;
770 }
771 m += 0x3FFF; /* biased exponent of 2^(M) */
772
773 expcont1:
774 fp2 = fp0; /* N */
775 fp0 = floatx80_mul(fp0, float32_to_floatx80(
776 make_float32(0xBC317218), status),
777 status); /* N * L1, L1 = lead(-log2/64) */
778 l2 = packFloatx80(0, 0x3FDC, UINT64_C(0x82E308654361C4C6));
779 fp2 = floatx80_mul(fp2, l2, status); /* N * L2, L1+L2 = -log2/64 */
780 fp0 = floatx80_add(fp0, fp1, status); /* X + N*L1 */
781 fp0 = floatx80_add(fp0, fp2, status); /* R */
782
783 fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */
784 fp2 = float32_to_floatx80(make_float32(0x3AB60B70),
785 status); /* A5 */
786 fp2 = floatx80_mul(fp2, fp1, status); /* fp2 is S*A5 */
787 fp3 = floatx80_mul(float32_to_floatx80(make_float32(0x3C088895),
788 status), fp1,
789 status); /* fp3 is S*A4 */
790 fp2 = floatx80_add(fp2, float64_to_floatx80(make_float64(
791 0x3FA5555555554431), status),
792 status); /* fp2 is A3+S*A5 */
793 fp3 = floatx80_add(fp3, float64_to_floatx80(make_float64(
794 0x3FC5555555554018), status),
795 status); /* fp3 is A2+S*A4 */
796 fp2 = floatx80_mul(fp2, fp1, status); /* fp2 is S*(A3+S*A5) */
797 fp3 = floatx80_mul(fp3, fp1, status); /* fp3 is S*(A2+S*A4) */
798 fp2 = floatx80_add(fp2, float32_to_floatx80(
799 make_float32(0x3F000000), status),
800 status); /* fp2 is A1+S*(A3+S*A5) */
801 fp3 = floatx80_mul(fp3, fp0, status); /* fp3 IS R*S*(A2+S*A4) */
802 fp2 = floatx80_mul(fp2, fp1,
803 status); /* fp2 IS S*(A1+S*(A3+S*A5)) */
804 fp0 = floatx80_add(fp0, fp3, status); /* fp0 IS R+R*S*(A2+S*A4) */
805 fp0 = floatx80_add(fp0, fp2, status); /* fp0 IS EXP(R) - 1 */
806
807 fp1 = exp_tbl[j];
808 fp0 = floatx80_mul(fp0, fp1, status); /* 2^(J/64)*(Exp(R)-1) */
809 fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j], status),
810 status); /* accurate 2^(J/64) */
811 fp0 = floatx80_add(fp0, fp1,
812 status); /* 2^(J/64) + 2^(J/64)*(Exp(R)-1) */
813
814 scale = packFloatx80(0, m, one_sig);
815 if (adjflag) {
816 adjscale = packFloatx80(0, m1, one_sig);
817 fp0 = floatx80_mul(fp0, adjscale, status);
818 }
819
820 status->float_rounding_mode = user_rnd_mode;
821 status->floatx80_rounding_precision = user_rnd_prec;
822
823 a = floatx80_mul(fp0, scale, status);
824
825 float_raise(float_flag_inexact, status);
826
827 return a;
828 } else { /* |X| >= 16380 log2 */
829 if (compact > 0x400CB27C) { /* |X| >= 16480 log2 */
830 status->float_rounding_mode = user_rnd_mode;
831 status->floatx80_rounding_precision = user_rnd_prec;
832 if (aSign) {
833 a = roundAndPackFloatx80(
834 status->floatx80_rounding_precision,
835 0, -0x1000, aSig, 0, status);
836 } else {
837 a = roundAndPackFloatx80(
838 status->floatx80_rounding_precision,
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 status->float_rounding_mode = user_rnd_mode;
880 status->floatx80_rounding_precision = user_rnd_prec;
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
floatx80_twotox(floatx80 a,float_status * status)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 = status->float_rounding_mode;
926 user_rnd_prec = status->floatx80_rounding_precision;
927 status->float_rounding_mode = float_round_nearest_even;
928 status->floatx80_rounding_precision = floatx80_precision_x;
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 status->float_rounding_mode = user_rnd_mode;
938 status->floatx80_rounding_precision = user_rnd_prec;
939
940 if (aSign) {
941 return roundAndPackFloatx80(status->floatx80_rounding_precision,
942 0, -0x1000, aSig, 0, status);
943 } else {
944 return roundAndPackFloatx80(status->floatx80_rounding_precision,
945 0, 0x8000, aSig, 0, status);
946 }
947 } else { /* |X| < 2^(-70) */
948 status->float_rounding_mode = user_rnd_mode;
949 status->floatx80_rounding_precision = user_rnd_prec;
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 status->float_rounding_mode = user_rnd_mode;
1031 status->floatx80_rounding_precision = user_rnd_prec;
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
floatx80_tentox(floatx80 a,float_status * status)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 = status->float_rounding_mode;
1076 user_rnd_prec = status->floatx80_rounding_precision;
1077 status->float_rounding_mode = float_round_nearest_even;
1078 status->floatx80_rounding_precision = floatx80_precision_x;
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 status->float_rounding_mode = user_rnd_mode;
1088 status->floatx80_rounding_precision = user_rnd_prec;
1089
1090 if (aSign) {
1091 return roundAndPackFloatx80(status->floatx80_rounding_precision,
1092 0, -0x1000, aSig, 0, status);
1093 } else {
1094 return roundAndPackFloatx80(status->floatx80_rounding_precision,
1095 0, 0x8000, aSig, 0, status);
1096 }
1097 } else { /* |X| < 2^(-70) */
1098 status->float_rounding_mode = user_rnd_mode;
1099 status->floatx80_rounding_precision = user_rnd_prec;
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 status->float_rounding_mode = user_rnd_mode;
1186 status->floatx80_rounding_precision = user_rnd_prec;
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
floatx80_tan(floatx80 a,float_status * status)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 = status->float_rounding_mode;
1231 user_rnd_prec = status->floatx80_rounding_precision;
1232 status->float_rounding_mode = float_round_nearest_even;
1233 status->floatx80_rounding_precision = floatx80_precision_x;
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 status->float_rounding_mode = user_rnd_mode;
1298 status->floatx80_rounding_precision = user_rnd_prec;
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 status->float_rounding_mode = user_rnd_mode;
1358 status->floatx80_rounding_precision = user_rnd_prec;
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 status->float_rounding_mode = user_rnd_mode;
1396 status->floatx80_rounding_precision = user_rnd_prec;
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
floatx80_sin(floatx80 a,float_status * status)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 = status->float_rounding_mode;
1442 user_rnd_prec = status->floatx80_rounding_precision;
1443 status->float_rounding_mode = float_round_nearest_even;
1444 status->floatx80_rounding_precision = floatx80_precision_x;
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 status->float_rounding_mode = user_rnd_mode;
1513 status->floatx80_rounding_precision = user_rnd_prec;
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 status->float_rounding_mode = user_rnd_mode;
1586 status->floatx80_rounding_precision = user_rnd_prec;
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 status->float_rounding_mode = user_rnd_mode;
1636 status->floatx80_rounding_precision = user_rnd_prec;
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
floatx80_cos(floatx80 a,float_status * status)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 = status->float_rounding_mode;
1682 user_rnd_prec = status->floatx80_rounding_precision;
1683 status->float_rounding_mode = float_round_nearest_even;
1684 status->floatx80_rounding_precision = floatx80_precision_x;
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 status->float_rounding_mode = user_rnd_mode;
1752 status->floatx80_rounding_precision = user_rnd_prec;
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 status->float_rounding_mode = user_rnd_mode;
1826 status->floatx80_rounding_precision = user_rnd_prec;
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 status->float_rounding_mode = user_rnd_mode;
1874 status->floatx80_rounding_precision = user_rnd_prec;
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
floatx80_atan(floatx80 a,float_status * status)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 = status->float_rounding_mode;
1921 user_rnd_prec = status->floatx80_rounding_precision;
1922 status->float_rounding_mode = float_round_nearest_even;
1923 status->floatx80_rounding_precision = floatx80_precision_x;
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 status->float_rounding_mode = user_rnd_mode;
1933 status->floatx80_rounding_precision = user_rnd_prec;
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 status->float_rounding_mode = user_rnd_mode;
1973 status->floatx80_rounding_precision = user_rnd_prec;
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 status->float_rounding_mode = user_rnd_mode;
1984 status->floatx80_rounding_precision = user_rnd_prec;
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 status->float_rounding_mode = user_rnd_mode;
2024 status->floatx80_rounding_precision = user_rnd_prec;
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 status->float_rounding_mode = user_rnd_mode;
2073 status->floatx80_rounding_precision = user_rnd_prec;
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
floatx80_asin(floatx80 a,float_status * status)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 = status->float_rounding_mode;
2126 user_rnd_prec = status->floatx80_rounding_precision;
2127 status->float_rounding_mode = float_round_nearest_even;
2128 status->floatx80_rounding_precision = floatx80_precision_x;
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 status->float_rounding_mode = user_rnd_mode;
2140 status->floatx80_rounding_precision = user_rnd_prec;
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
floatx80_acos(floatx80 a,float_status * status)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(status->floatx80_rounding_precision, 0,
2175 piby2_exp, pi_sig, 0, status);
2176 }
2177
2178 compact = floatx80_make_compact(aExp, aSig);
2179
2180 if (compact >= 0x3FFF8000) { /* |X| >= 1 */
2181 if (aExp == one_exp && aSig == one_sig) { /* |X| == 1 */
2182 if (aSign) { /* X == -1 */
2183 a = packFloatx80(0, pi_exp, pi_sig);
2184 float_raise(float_flag_inexact, status);
2185 return floatx80_move(a, status);
2186 } else { /* X == +1 */
2187 return packFloatx80(0, 0, 0);
2188 }
2189 } else { /* |X| > 1 */
2190 float_raise(float_flag_invalid, status);
2191 return floatx80_default_nan(status);
2192 }
2193 } /* |X| < 1 */
2194
2195 user_rnd_mode = status->float_rounding_mode;
2196 user_rnd_prec = status->floatx80_rounding_precision;
2197 status->float_rounding_mode = float_round_nearest_even;
2198 status->floatx80_rounding_precision = floatx80_precision_x;
2199
2200 one = packFloatx80(0, one_exp, one_sig);
2201 fp0 = a;
2202
2203 fp1 = floatx80_add(one, fp0, status); /* 1 + X */
2204 fp0 = floatx80_sub(one, fp0, status); /* 1 - X */
2205 fp0 = floatx80_div(fp0, fp1, status); /* (1-X)/(1+X) */
2206 fp0 = floatx80_sqrt(fp0, status); /* SQRT((1-X)/(1+X)) */
2207 fp0 = floatx80_atan(fp0, status); /* ATAN(SQRT((1-X)/(1+X))) */
2208
2209 status->float_rounding_mode = user_rnd_mode;
2210 status->floatx80_rounding_precision = user_rnd_prec;
2211
2212 a = floatx80_add(fp0, fp0, status); /* 2 * ATAN(SQRT((1-X)/(1+X))) */
2213
2214 float_raise(float_flag_inexact, status);
2215
2216 return a;
2217 }
2218
2219 /*
2220 * Hyperbolic arc tangent
2221 */
2222
floatx80_atanh(floatx80 a,float_status * status)2223 floatx80 floatx80_atanh(floatx80 a, float_status *status)
2224 {
2225 bool aSign;
2226 int32_t aExp;
2227 uint64_t aSig;
2228
2229 FloatRoundMode user_rnd_mode;
2230 FloatX80RoundPrec user_rnd_prec;
2231
2232 int32_t compact;
2233 floatx80 fp0, fp1, fp2, one;
2234
2235 aSig = extractFloatx80Frac(a);
2236 aExp = extractFloatx80Exp(a);
2237 aSign = extractFloatx80Sign(a);
2238
2239 if (aExp == 0x7FFF && (uint64_t) (aSig << 1)) {
2240 return propagateFloatx80NaNOneArg(a, status);
2241 }
2242
2243 if (aExp == 0 && aSig == 0) {
2244 return packFloatx80(aSign, 0, 0);
2245 }
2246
2247 compact = floatx80_make_compact(aExp, aSig);
2248
2249 if (compact >= 0x3FFF8000) { /* |X| >= 1 */
2250 if (aExp == one_exp && aSig == one_sig) { /* |X| == 1 */
2251 float_raise(float_flag_divbyzero, status);
2252 return floatx80_default_inf(aSign, status);
2253 } else { /* |X| > 1 */
2254 float_raise(float_flag_invalid, status);
2255 return floatx80_default_nan(status);
2256 }
2257 } /* |X| < 1 */
2258
2259 user_rnd_mode = status->float_rounding_mode;
2260 user_rnd_prec = status->floatx80_rounding_precision;
2261 status->float_rounding_mode = float_round_nearest_even;
2262 status->floatx80_rounding_precision = floatx80_precision_x;
2263
2264 one = packFloatx80(0, one_exp, one_sig);
2265 fp2 = packFloatx80(aSign, 0x3FFE, one_sig); /* SIGN(X) * (1/2) */
2266 fp0 = packFloatx80(0, aExp, aSig); /* Y = |X| */
2267 fp1 = packFloatx80(1, aExp, aSig); /* -Y */
2268 fp0 = floatx80_add(fp0, fp0, status); /* 2Y */
2269 fp1 = floatx80_add(fp1, one, status); /* 1-Y */
2270 fp0 = floatx80_div(fp0, fp1, status); /* Z = 2Y/(1-Y) */
2271 fp0 = floatx80_lognp1(fp0, status); /* LOG1P(Z) */
2272
2273 status->float_rounding_mode = user_rnd_mode;
2274 status->floatx80_rounding_precision = user_rnd_prec;
2275
2276 a = floatx80_mul(fp0, fp2,
2277 status); /* ATANH(X) = SIGN(X) * (1/2) * LOG1P(Z) */
2278
2279 float_raise(float_flag_inexact, status);
2280
2281 return a;
2282 }
2283
2284 /*
2285 * e to x minus 1
2286 */
2287
floatx80_etoxm1(floatx80 a,float_status * status)2288 floatx80 floatx80_etoxm1(floatx80 a, float_status *status)
2289 {
2290 bool aSign;
2291 int32_t aExp;
2292 uint64_t aSig;
2293
2294 FloatRoundMode user_rnd_mode;
2295 FloatX80RoundPrec user_rnd_prec;
2296
2297 int32_t compact, n, j, m, m1;
2298 floatx80 fp0, fp1, fp2, fp3, l2, sc, onebysc;
2299
2300 aSig = extractFloatx80Frac(a);
2301 aExp = extractFloatx80Exp(a);
2302 aSign = extractFloatx80Sign(a);
2303
2304 if (aExp == 0x7FFF) {
2305 if ((uint64_t) (aSig << 1)) {
2306 return propagateFloatx80NaNOneArg(a, status);
2307 }
2308 if (aSign) {
2309 return packFloatx80(aSign, one_exp, one_sig);
2310 }
2311 return floatx80_default_inf(0, status);
2312 }
2313
2314 if (aExp == 0 && aSig == 0) {
2315 return packFloatx80(aSign, 0, 0);
2316 }
2317
2318 user_rnd_mode = status->float_rounding_mode;
2319 user_rnd_prec = status->floatx80_rounding_precision;
2320 status->float_rounding_mode = float_round_nearest_even;
2321 status->floatx80_rounding_precision = floatx80_precision_x;
2322
2323 if (aExp >= 0x3FFD) { /* |X| >= 1/4 */
2324 compact = floatx80_make_compact(aExp, aSig);
2325
2326 if (compact <= 0x4004C215) { /* |X| <= 70 log2 */
2327 fp0 = a;
2328 fp1 = a;
2329 fp0 = floatx80_mul(fp0, float32_to_floatx80(
2330 make_float32(0x42B8AA3B), status),
2331 status); /* 64/log2 * X */
2332 n = floatx80_to_int32(fp0, status); /* int(64/log2*X) */
2333 fp0 = int32_to_floatx80(n, status);
2334
2335 j = n & 0x3F; /* J = N mod 64 */
2336 m = n / 64; /* NOTE: this is really arithmetic right shift by 6 */
2337 if (n < 0 && j) {
2338 /*
2339 * arithmetic right shift is division and
2340 * round towards minus infinity
2341 */
2342 m--;
2343 }
2344 m1 = -m;
2345 /*m += 0x3FFF; // biased exponent of 2^(M) */
2346 /*m1 += 0x3FFF; // biased exponent of -2^(-M) */
2347
2348 fp2 = fp0; /* N */
2349 fp0 = floatx80_mul(fp0, float32_to_floatx80(
2350 make_float32(0xBC317218), status),
2351 status); /* N * L1, L1 = lead(-log2/64) */
2352 l2 = packFloatx80(0, 0x3FDC, UINT64_C(0x82E308654361C4C6));
2353 fp2 = floatx80_mul(fp2, l2, status); /* N * L2, L1+L2 = -log2/64 */
2354 fp0 = floatx80_add(fp0, fp1, status); /* X + N*L1 */
2355 fp0 = floatx80_add(fp0, fp2, status); /* R */
2356
2357 fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */
2358 fp2 = float32_to_floatx80(make_float32(0x3950097B),
2359 status); /* A6 */
2360 fp2 = floatx80_mul(fp2, fp1, status); /* fp2 is S*A6 */
2361 fp3 = floatx80_mul(float32_to_floatx80(make_float32(0x3AB60B6A),
2362 status), fp1, status); /* fp3 is S*A5 */
2363 fp2 = floatx80_add(fp2, float64_to_floatx80(
2364 make_float64(0x3F81111111174385), status),
2365 status); /* fp2 IS A4+S*A6 */
2366 fp3 = floatx80_add(fp3, float64_to_floatx80(
2367 make_float64(0x3FA5555555554F5A), status),
2368 status); /* fp3 is A3+S*A5 */
2369 fp2 = floatx80_mul(fp2, fp1, status); /* fp2 IS S*(A4+S*A6) */
2370 fp3 = floatx80_mul(fp3, fp1, status); /* fp3 IS S*(A3+S*A5) */
2371 fp2 = floatx80_add(fp2, float64_to_floatx80(
2372 make_float64(0x3FC5555555555555), status),
2373 status); /* fp2 IS A2+S*(A4+S*A6) */
2374 fp3 = floatx80_add(fp3, float32_to_floatx80(
2375 make_float32(0x3F000000), status),
2376 status); /* fp3 IS A1+S*(A3+S*A5) */
2377 fp2 = floatx80_mul(fp2, fp1,
2378 status); /* fp2 IS S*(A2+S*(A4+S*A6)) */
2379 fp1 = floatx80_mul(fp1, fp3,
2380 status); /* fp1 IS S*(A1+S*(A3+S*A5)) */
2381 fp2 = floatx80_mul(fp2, fp0,
2382 status); /* fp2 IS R*S*(A2+S*(A4+S*A6)) */
2383 fp0 = floatx80_add(fp0, fp1,
2384 status); /* fp0 IS R+S*(A1+S*(A3+S*A5)) */
2385 fp0 = floatx80_add(fp0, fp2, status); /* fp0 IS EXP(R) - 1 */
2386
2387 fp0 = floatx80_mul(fp0, exp_tbl[j],
2388 status); /* 2^(J/64)*(Exp(R)-1) */
2389
2390 if (m >= 64) {
2391 fp1 = float32_to_floatx80(exp_tbl2[j], status);
2392 onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); /* -2^(-M) */
2393 fp1 = floatx80_add(fp1, onebysc, status);
2394 fp0 = floatx80_add(fp0, fp1, status);
2395 fp0 = floatx80_add(fp0, exp_tbl[j], status);
2396 } else if (m < -3) {
2397 fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j],
2398 status), status);
2399 fp0 = floatx80_add(fp0, exp_tbl[j], status);
2400 onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); /* -2^(-M) */
2401 fp0 = floatx80_add(fp0, onebysc, status);
2402 } else { /* -3 <= m <= 63 */
2403 fp1 = exp_tbl[j];
2404 fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j],
2405 status), status);
2406 onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); /* -2^(-M) */
2407 fp1 = floatx80_add(fp1, onebysc, status);
2408 fp0 = floatx80_add(fp0, fp1, status);
2409 }
2410
2411 sc = packFloatx80(0, m + 0x3FFF, one_sig);
2412
2413 status->float_rounding_mode = user_rnd_mode;
2414 status->floatx80_rounding_precision = user_rnd_prec;
2415
2416 a = floatx80_mul(fp0, sc, status);
2417
2418 float_raise(float_flag_inexact, status);
2419
2420 return a;
2421 } else { /* |X| > 70 log2 */
2422 if (aSign) {
2423 fp0 = float32_to_floatx80(make_float32(0xBF800000),
2424 status); /* -1 */
2425
2426 status->float_rounding_mode = user_rnd_mode;
2427 status->floatx80_rounding_precision = user_rnd_prec;
2428
2429 a = floatx80_add(fp0, float32_to_floatx80(
2430 make_float32(0x00800000), status),
2431 status); /* -1 + 2^(-126) */
2432
2433 float_raise(float_flag_inexact, status);
2434
2435 return a;
2436 } else {
2437 status->float_rounding_mode = user_rnd_mode;
2438 status->floatx80_rounding_precision = user_rnd_prec;
2439
2440 return floatx80_etox(a, status);
2441 }
2442 }
2443 } else { /* |X| < 1/4 */
2444 if (aExp >= 0x3FBE) {
2445 fp0 = a;
2446 fp0 = floatx80_mul(fp0, fp0, status); /* S = X*X */
2447 fp1 = float32_to_floatx80(make_float32(0x2F30CAA8),
2448 status); /* B12 */
2449 fp1 = floatx80_mul(fp1, fp0, status); /* S * B12 */
2450 fp2 = float32_to_floatx80(make_float32(0x310F8290),
2451 status); /* B11 */
2452 fp1 = floatx80_add(fp1, float32_to_floatx80(
2453 make_float32(0x32D73220), status),
2454 status); /* B10 */
2455 fp2 = floatx80_mul(fp2, fp0, status);
2456 fp1 = floatx80_mul(fp1, fp0, status);
2457 fp2 = floatx80_add(fp2, float32_to_floatx80(
2458 make_float32(0x3493F281), status),
2459 status); /* B9 */
2460 fp1 = floatx80_add(fp1, float64_to_floatx80(
2461 make_float64(0x3EC71DE3A5774682), status),
2462 status); /* B8 */
2463 fp2 = floatx80_mul(fp2, fp0, status);
2464 fp1 = floatx80_mul(fp1, fp0, status);
2465 fp2 = floatx80_add(fp2, float64_to_floatx80(
2466 make_float64(0x3EFA01A019D7CB68), status),
2467 status); /* B7 */
2468 fp1 = floatx80_add(fp1, float64_to_floatx80(
2469 make_float64(0x3F2A01A01A019DF3), status),
2470 status); /* B6 */
2471 fp2 = floatx80_mul(fp2, fp0, status);
2472 fp1 = floatx80_mul(fp1, fp0, status);
2473 fp2 = floatx80_add(fp2, float64_to_floatx80(
2474 make_float64(0x3F56C16C16C170E2), status),
2475 status); /* B5 */
2476 fp1 = floatx80_add(fp1, float64_to_floatx80(
2477 make_float64(0x3F81111111111111), status),
2478 status); /* B4 */
2479 fp2 = floatx80_mul(fp2, fp0, status);
2480 fp1 = floatx80_mul(fp1, fp0, status);
2481 fp2 = floatx80_add(fp2, float64_to_floatx80(
2482 make_float64(0x3FA5555555555555), status),
2483 status); /* B3 */
2484 fp3 = packFloatx80(0, 0x3FFC, UINT64_C(0xAAAAAAAAAAAAAAAB));
2485 fp1 = floatx80_add(fp1, fp3, status); /* B2 */
2486 fp2 = floatx80_mul(fp2, fp0, status);
2487 fp1 = floatx80_mul(fp1, fp0, status);
2488
2489 fp2 = floatx80_mul(fp2, fp0, status);
2490 fp1 = floatx80_mul(fp1, a, status);
2491
2492 fp0 = floatx80_mul(fp0, float32_to_floatx80(
2493 make_float32(0x3F000000), status),
2494 status); /* S*B1 */
2495 fp1 = floatx80_add(fp1, fp2, status); /* Q */
2496 fp0 = floatx80_add(fp0, fp1, status); /* S*B1+Q */
2497
2498 status->float_rounding_mode = user_rnd_mode;
2499 status->floatx80_rounding_precision = user_rnd_prec;
2500
2501 a = floatx80_add(fp0, a, status);
2502
2503 float_raise(float_flag_inexact, status);
2504
2505 return a;
2506 } else { /* |X| < 2^(-65) */
2507 sc = packFloatx80(1, 1, one_sig);
2508 fp0 = a;
2509
2510 if (aExp < 0x0033) { /* |X| < 2^(-16382) */
2511 fp0 = floatx80_mul(fp0, float64_to_floatx80(
2512 make_float64(0x48B0000000000000), status),
2513 status);
2514 fp0 = floatx80_add(fp0, sc, status);
2515
2516 status->float_rounding_mode = user_rnd_mode;
2517 status->floatx80_rounding_precision = user_rnd_prec;
2518
2519 a = floatx80_mul(fp0, float64_to_floatx80(
2520 make_float64(0x3730000000000000), status),
2521 status);
2522 } else {
2523 status->float_rounding_mode = user_rnd_mode;
2524 status->floatx80_rounding_precision = user_rnd_prec;
2525
2526 a = floatx80_add(fp0, sc, status);
2527 }
2528
2529 float_raise(float_flag_inexact, status);
2530
2531 return a;
2532 }
2533 }
2534 }
2535
2536 /*
2537 * Hyperbolic tangent
2538 */
2539
floatx80_tanh(floatx80 a,float_status * status)2540 floatx80 floatx80_tanh(floatx80 a, float_status *status)
2541 {
2542 bool aSign, vSign;
2543 int32_t aExp, vExp;
2544 uint64_t aSig, vSig;
2545
2546 FloatRoundMode user_rnd_mode;
2547 FloatX80RoundPrec user_rnd_prec;
2548
2549 int32_t compact;
2550 floatx80 fp0, fp1;
2551 uint32_t sign;
2552
2553 aSig = extractFloatx80Frac(a);
2554 aExp = extractFloatx80Exp(a);
2555 aSign = extractFloatx80Sign(a);
2556
2557 if (aExp == 0x7FFF) {
2558 if ((uint64_t) (aSig << 1)) {
2559 return propagateFloatx80NaNOneArg(a, status);
2560 }
2561 return packFloatx80(aSign, one_exp, one_sig);
2562 }
2563
2564 if (aExp == 0 && aSig == 0) {
2565 return packFloatx80(aSign, 0, 0);
2566 }
2567
2568 user_rnd_mode = status->float_rounding_mode;
2569 user_rnd_prec = status->floatx80_rounding_precision;
2570 status->float_rounding_mode = float_round_nearest_even;
2571 status->floatx80_rounding_precision = floatx80_precision_x;
2572
2573 compact = floatx80_make_compact(aExp, aSig);
2574
2575 if (compact < 0x3FD78000 || compact > 0x3FFFDDCE) {
2576 /* TANHBORS */
2577 if (compact < 0x3FFF8000) {
2578 /* TANHSM */
2579 status->float_rounding_mode = user_rnd_mode;
2580 status->floatx80_rounding_precision = user_rnd_prec;
2581
2582 a = floatx80_move(a, status);
2583
2584 float_raise(float_flag_inexact, status);
2585
2586 return a;
2587 } else {
2588 if (compact > 0x40048AA1) {
2589 /* TANHHUGE */
2590 sign = 0x3F800000;
2591 sign |= aSign ? 0x80000000 : 0x00000000;
2592 fp0 = float32_to_floatx80(make_float32(sign), status);
2593 sign &= 0x80000000;
2594 sign ^= 0x80800000; /* -SIGN(X)*EPS */
2595
2596 status->float_rounding_mode = user_rnd_mode;
2597 status->floatx80_rounding_precision = user_rnd_prec;
2598
2599 a = floatx80_add(fp0, float32_to_floatx80(make_float32(sign),
2600 status), status);
2601
2602 float_raise(float_flag_inexact, status);
2603
2604 return a;
2605 } else {
2606 fp0 = packFloatx80(0, aExp + 1, aSig); /* Y = 2|X| */
2607 fp0 = floatx80_etox(fp0, status); /* FP0 IS EXP(Y) */
2608 fp0 = floatx80_add(fp0, float32_to_floatx80(
2609 make_float32(0x3F800000),
2610 status), status); /* EXP(Y)+1 */
2611 sign = aSign ? 0x80000000 : 0x00000000;
2612 fp1 = floatx80_div(float32_to_floatx80(make_float32(
2613 sign ^ 0xC0000000), status), fp0,
2614 status); /* -SIGN(X)*2 / [EXP(Y)+1] */
2615 fp0 = float32_to_floatx80(make_float32(sign | 0x3F800000),
2616 status); /* SIGN */
2617
2618 status->float_rounding_mode = user_rnd_mode;
2619 status->floatx80_rounding_precision = user_rnd_prec;
2620
2621 a = floatx80_add(fp1, fp0, status);
2622
2623 float_raise(float_flag_inexact, status);
2624
2625 return a;
2626 }
2627 }
2628 } else { /* 2**(-40) < |X| < (5/2)LOG2 */
2629 fp0 = packFloatx80(0, aExp + 1, aSig); /* Y = 2|X| */
2630 fp0 = floatx80_etoxm1(fp0, status); /* FP0 IS Z = EXPM1(Y) */
2631 fp1 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x40000000),
2632 status),
2633 status); /* Z+2 */
2634
2635 vSign = extractFloatx80Sign(fp1);
2636 vExp = extractFloatx80Exp(fp1);
2637 vSig = extractFloatx80Frac(fp1);
2638
2639 fp1 = packFloatx80(vSign ^ aSign, vExp, vSig);
2640
2641 status->float_rounding_mode = user_rnd_mode;
2642 status->floatx80_rounding_precision = user_rnd_prec;
2643
2644 a = floatx80_div(fp0, fp1, status);
2645
2646 float_raise(float_flag_inexact, status);
2647
2648 return a;
2649 }
2650 }
2651
2652 /*
2653 * Hyperbolic sine
2654 */
2655
floatx80_sinh(floatx80 a,float_status * status)2656 floatx80 floatx80_sinh(floatx80 a, float_status *status)
2657 {
2658 bool aSign;
2659 int32_t aExp;
2660 uint64_t aSig;
2661
2662 FloatRoundMode user_rnd_mode;
2663 FloatX80RoundPrec user_rnd_prec;
2664
2665 int32_t compact;
2666 floatx80 fp0, fp1, fp2;
2667 float32 fact;
2668
2669 aSig = extractFloatx80Frac(a);
2670 aExp = extractFloatx80Exp(a);
2671 aSign = extractFloatx80Sign(a);
2672
2673 if (aExp == 0x7FFF) {
2674 if ((uint64_t) (aSig << 1)) {
2675 return propagateFloatx80NaNOneArg(a, status);
2676 }
2677 return floatx80_default_inf(aSign, status);
2678 }
2679
2680 if (aExp == 0 && aSig == 0) {
2681 return packFloatx80(aSign, 0, 0);
2682 }
2683
2684 user_rnd_mode = status->float_rounding_mode;
2685 user_rnd_prec = status->floatx80_rounding_precision;
2686 status->float_rounding_mode = float_round_nearest_even;
2687 status->floatx80_rounding_precision = floatx80_precision_x;
2688
2689 compact = floatx80_make_compact(aExp, aSig);
2690
2691 if (compact > 0x400CB167) {
2692 /* SINHBIG */
2693 if (compact > 0x400CB2B3) {
2694 status->float_rounding_mode = user_rnd_mode;
2695 status->floatx80_rounding_precision = user_rnd_prec;
2696
2697 return roundAndPackFloatx80(status->floatx80_rounding_precision,
2698 aSign, 0x8000, aSig, 0, status);
2699 } else {
2700 fp0 = floatx80_abs(a); /* Y = |X| */
2701 fp0 = floatx80_sub(fp0, float64_to_floatx80(
2702 make_float64(0x40C62D38D3D64634), status),
2703 status); /* (|X|-16381LOG2_LEAD) */
2704 fp0 = floatx80_sub(fp0, float64_to_floatx80(
2705 make_float64(0x3D6F90AEB1E75CC7), status),
2706 status); /* |X| - 16381 LOG2, ACCURATE */
2707 fp0 = floatx80_etox(fp0, status);
2708 fp2 = packFloatx80(aSign, 0x7FFB, one_sig);
2709
2710 status->float_rounding_mode = user_rnd_mode;
2711 status->floatx80_rounding_precision = user_rnd_prec;
2712
2713 a = floatx80_mul(fp0, fp2, status);
2714
2715 float_raise(float_flag_inexact, status);
2716
2717 return a;
2718 }
2719 } else { /* |X| < 16380 LOG2 */
2720 fp0 = floatx80_abs(a); /* Y = |X| */
2721 fp0 = floatx80_etoxm1(fp0, status); /* FP0 IS Z = EXPM1(Y) */
2722 fp1 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000),
2723 status), status); /* 1+Z */
2724 fp2 = fp0;
2725 fp0 = floatx80_div(fp0, fp1, status); /* Z/(1+Z) */
2726 fp0 = floatx80_add(fp0, fp2, status);
2727
2728 fact = packFloat32(aSign, 0x7E, 0);
2729
2730 status->float_rounding_mode = user_rnd_mode;
2731 status->floatx80_rounding_precision = user_rnd_prec;
2732
2733 a = floatx80_mul(fp0, float32_to_floatx80(fact, status), status);
2734
2735 float_raise(float_flag_inexact, status);
2736
2737 return a;
2738 }
2739 }
2740
2741 /*
2742 * Hyperbolic cosine
2743 */
2744
floatx80_cosh(floatx80 a,float_status * status)2745 floatx80 floatx80_cosh(floatx80 a, float_status *status)
2746 {
2747 int32_t aExp;
2748 uint64_t aSig;
2749
2750 FloatRoundMode user_rnd_mode;
2751 FloatX80RoundPrec user_rnd_prec;
2752
2753 int32_t compact;
2754 floatx80 fp0, fp1;
2755
2756 aSig = extractFloatx80Frac(a);
2757 aExp = extractFloatx80Exp(a);
2758
2759 if (aExp == 0x7FFF) {
2760 if ((uint64_t) (aSig << 1)) {
2761 return propagateFloatx80NaNOneArg(a, status);
2762 }
2763 return floatx80_default_inf(0, status);
2764 }
2765
2766 if (aExp == 0 && aSig == 0) {
2767 return packFloatx80(0, one_exp, one_sig);
2768 }
2769
2770 user_rnd_mode = status->float_rounding_mode;
2771 user_rnd_prec = status->floatx80_rounding_precision;
2772 status->float_rounding_mode = float_round_nearest_even;
2773 status->floatx80_rounding_precision = floatx80_precision_x;
2774
2775 compact = floatx80_make_compact(aExp, aSig);
2776
2777 if (compact > 0x400CB167) {
2778 if (compact > 0x400CB2B3) {
2779 status->float_rounding_mode = user_rnd_mode;
2780 status->floatx80_rounding_precision = user_rnd_prec;
2781 return roundAndPackFloatx80(status->floatx80_rounding_precision, 0,
2782 0x8000, one_sig, 0, status);
2783 } else {
2784 fp0 = packFloatx80(0, aExp, aSig);
2785 fp0 = floatx80_sub(fp0, float64_to_floatx80(
2786 make_float64(0x40C62D38D3D64634), status),
2787 status);
2788 fp0 = floatx80_sub(fp0, float64_to_floatx80(
2789 make_float64(0x3D6F90AEB1E75CC7), status),
2790 status);
2791 fp0 = floatx80_etox(fp0, status);
2792 fp1 = packFloatx80(0, 0x7FFB, one_sig);
2793
2794 status->float_rounding_mode = user_rnd_mode;
2795 status->floatx80_rounding_precision = user_rnd_prec;
2796
2797 a = floatx80_mul(fp0, fp1, status);
2798
2799 float_raise(float_flag_inexact, status);
2800
2801 return a;
2802 }
2803 }
2804
2805 fp0 = packFloatx80(0, aExp, aSig); /* |X| */
2806 fp0 = floatx80_etox(fp0, status); /* EXP(|X|) */
2807 fp0 = floatx80_mul(fp0, float32_to_floatx80(make_float32(0x3F000000),
2808 status), status); /* (1/2)*EXP(|X|) */
2809 fp1 = float32_to_floatx80(make_float32(0x3E800000), status); /* 1/4 */
2810 fp1 = floatx80_div(fp1, fp0, status); /* 1/(2*EXP(|X|)) */
2811
2812 status->float_rounding_mode = user_rnd_mode;
2813 status->floatx80_rounding_precision = user_rnd_prec;
2814
2815 a = floatx80_add(fp0, fp1, status);
2816
2817 float_raise(float_flag_inexact, status);
2818
2819 return a;
2820 }
2821