Loading...
1// SPDX-License-Identifier: GPL-2.0
2/*---------------------------------------------------------------------------+
3 | fpu_trig.c |
4 | |
5 | Implementation of the FPU "transcendental" functions. |
6 | |
7 | Copyright (C) 1992,1993,1994,1997,1999 |
8 | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, |
9 | Australia. E-mail billm@melbpc.org.au |
10 | |
11 | |
12 +---------------------------------------------------------------------------*/
13
14#include "fpu_system.h"
15#include "exception.h"
16#include "fpu_emu.h"
17#include "status_w.h"
18#include "control_w.h"
19#include "reg_constant.h"
20
21static void rem_kernel(unsigned long long st0, unsigned long long *y,
22 unsigned long long st1, unsigned long long q, int n);
23
24#define BETTER_THAN_486
25
26#define FCOS 4
27
28/* Used only by fptan, fsin, fcos, and fsincos. */
29/* This routine produces very accurate results, similar to
30 using a value of pi with more than 128 bits precision. */
31/* Limited measurements show no results worse than 64 bit precision
32 except for the results for arguments close to 2^63, where the
33 precision of the result sometimes degrades to about 63.9 bits */
34static int trig_arg(FPU_REG *st0_ptr, int even)
35{
36 FPU_REG tmp;
37 u_char tmptag;
38 unsigned long long q;
39 int old_cw = control_word, saved_status = partial_status;
40 int tag, st0_tag = TAG_Valid;
41
42 if (exponent(st0_ptr) >= 63) {
43 partial_status |= SW_C2; /* Reduction incomplete. */
44 return -1;
45 }
46
47 control_word &= ~CW_RC;
48 control_word |= RC_CHOP;
49
50 setpositive(st0_ptr);
51 tag = FPU_u_div(st0_ptr, &CONST_PI2, &tmp, PR_64_BITS | RC_CHOP | 0x3f,
52 SIGN_POS);
53
54 FPU_round_to_int(&tmp, tag); /* Fortunately, this can't overflow
55 to 2^64 */
56 q = significand(&tmp);
57 if (q) {
58 rem_kernel(significand(st0_ptr),
59 &significand(&tmp),
60 significand(&CONST_PI2),
61 q, exponent(st0_ptr) - exponent(&CONST_PI2));
62 setexponent16(&tmp, exponent(&CONST_PI2));
63 st0_tag = FPU_normalize(&tmp);
64 FPU_copy_to_reg0(&tmp, st0_tag);
65 }
66
67 if ((even && !(q & 1)) || (!even && (q & 1))) {
68 st0_tag =
69 FPU_sub(REV | LOADED | TAG_Valid, (int)&CONST_PI2,
70 FULL_PRECISION);
71
72#ifdef BETTER_THAN_486
73 /* So far, the results are exact but based upon a 64 bit
74 precision approximation to pi/2. The technique used
75 now is equivalent to using an approximation to pi/2 which
76 is accurate to about 128 bits. */
77 if ((exponent(st0_ptr) <= exponent(&CONST_PI2extra) + 64)
78 || (q > 1)) {
79 /* This code gives the effect of having pi/2 to better than
80 128 bits precision. */
81
82 significand(&tmp) = q + 1;
83 setexponent16(&tmp, 63);
84 FPU_normalize(&tmp);
85 tmptag =
86 FPU_u_mul(&CONST_PI2extra, &tmp, &tmp,
87 FULL_PRECISION, SIGN_POS,
88 exponent(&CONST_PI2extra) +
89 exponent(&tmp));
90 setsign(&tmp, getsign(&CONST_PI2extra));
91 st0_tag = FPU_add(&tmp, tmptag, 0, FULL_PRECISION);
92 if (signnegative(st0_ptr)) {
93 /* CONST_PI2extra is negative, so the result of the addition
94 can be negative. This means that the argument is actually
95 in a different quadrant. The correction is always < pi/2,
96 so it can't overflow into yet another quadrant. */
97 setpositive(st0_ptr);
98 q++;
99 }
100 }
101#endif /* BETTER_THAN_486 */
102 }
103#ifdef BETTER_THAN_486
104 else {
105 /* So far, the results are exact but based upon a 64 bit
106 precision approximation to pi/2. The technique used
107 now is equivalent to using an approximation to pi/2 which
108 is accurate to about 128 bits. */
109 if (((q > 0)
110 && (exponent(st0_ptr) <= exponent(&CONST_PI2extra) + 64))
111 || (q > 1)) {
112 /* This code gives the effect of having p/2 to better than
113 128 bits precision. */
114
115 significand(&tmp) = q;
116 setexponent16(&tmp, 63);
117 FPU_normalize(&tmp); /* This must return TAG_Valid */
118 tmptag =
119 FPU_u_mul(&CONST_PI2extra, &tmp, &tmp,
120 FULL_PRECISION, SIGN_POS,
121 exponent(&CONST_PI2extra) +
122 exponent(&tmp));
123 setsign(&tmp, getsign(&CONST_PI2extra));
124 st0_tag = FPU_sub(LOADED | (tmptag & 0x0f), (int)&tmp,
125 FULL_PRECISION);
126 if ((exponent(st0_ptr) == exponent(&CONST_PI2)) &&
127 ((st0_ptr->sigh > CONST_PI2.sigh)
128 || ((st0_ptr->sigh == CONST_PI2.sigh)
129 && (st0_ptr->sigl > CONST_PI2.sigl)))) {
130 /* CONST_PI2extra is negative, so the result of the
131 subtraction can be larger than pi/2. This means
132 that the argument is actually in a different quadrant.
133 The correction is always < pi/2, so it can't overflow
134 into yet another quadrant. */
135 st0_tag =
136 FPU_sub(REV | LOADED | TAG_Valid,
137 (int)&CONST_PI2, FULL_PRECISION);
138 q++;
139 }
140 }
141 }
142#endif /* BETTER_THAN_486 */
143
144 FPU_settag0(st0_tag);
145 control_word = old_cw;
146 partial_status = saved_status & ~SW_C2; /* Reduction complete. */
147
148 return (q & 3) | even;
149}
150
151/* Convert a long to register */
152static void convert_l2reg(long const *arg, int deststnr)
153{
154 int tag;
155 long num = *arg;
156 u_char sign;
157 FPU_REG *dest = &st(deststnr);
158
159 if (num == 0) {
160 FPU_copy_to_regi(&CONST_Z, TAG_Zero, deststnr);
161 return;
162 }
163
164 if (num > 0) {
165 sign = SIGN_POS;
166 } else {
167 num = -num;
168 sign = SIGN_NEG;
169 }
170
171 dest->sigh = num;
172 dest->sigl = 0;
173 setexponent16(dest, 31);
174 tag = FPU_normalize(dest);
175 FPU_settagi(deststnr, tag);
176 setsign(dest, sign);
177 return;
178}
179
180static void single_arg_error(FPU_REG *st0_ptr, u_char st0_tag)
181{
182 if (st0_tag == TAG_Empty)
183 FPU_stack_underflow(); /* Puts a QNaN in st(0) */
184 else if (st0_tag == TW_NaN)
185 real_1op_NaN(st0_ptr); /* return with a NaN in st(0) */
186#ifdef PARANOID
187 else
188 EXCEPTION(EX_INTERNAL | 0x0112);
189#endif /* PARANOID */
190}
191
192static void single_arg_2_error(FPU_REG *st0_ptr, u_char st0_tag)
193{
194 int isNaN;
195
196 switch (st0_tag) {
197 case TW_NaN:
198 isNaN = (exponent(st0_ptr) == EXP_OVER)
199 && (st0_ptr->sigh & 0x80000000);
200 if (isNaN && !(st0_ptr->sigh & 0x40000000)) { /* Signaling ? */
201 EXCEPTION(EX_Invalid);
202 if (control_word & CW_Invalid) {
203 /* The masked response */
204 /* Convert to a QNaN */
205 st0_ptr->sigh |= 0x40000000;
206 push();
207 FPU_copy_to_reg0(st0_ptr, TAG_Special);
208 }
209 } else if (isNaN) {
210 /* A QNaN */
211 push();
212 FPU_copy_to_reg0(st0_ptr, TAG_Special);
213 } else {
214 /* pseudoNaN or other unsupported */
215 EXCEPTION(EX_Invalid);
216 if (control_word & CW_Invalid) {
217 /* The masked response */
218 FPU_copy_to_reg0(&CONST_QNaN, TAG_Special);
219 push();
220 FPU_copy_to_reg0(&CONST_QNaN, TAG_Special);
221 }
222 }
223 break; /* return with a NaN in st(0) */
224#ifdef PARANOID
225 default:
226 EXCEPTION(EX_INTERNAL | 0x0112);
227#endif /* PARANOID */
228 }
229}
230
231/*---------------------------------------------------------------------------*/
232
233static void f2xm1(FPU_REG *st0_ptr, u_char tag)
234{
235 FPU_REG a;
236
237 clear_C1();
238
239 if (tag == TAG_Valid) {
240 /* For an 80486 FPU, the result is undefined if the arg is >= 1.0 */
241 if (exponent(st0_ptr) < 0) {
242 denormal_arg:
243
244 FPU_to_exp16(st0_ptr, &a);
245
246 /* poly_2xm1(x) requires 0 < st(0) < 1. */
247 poly_2xm1(getsign(st0_ptr), &a, st0_ptr);
248 }
249 set_precision_flag_up(); /* 80486 appears to always do this */
250 return;
251 }
252
253 if (tag == TAG_Zero)
254 return;
255
256 if (tag == TAG_Special)
257 tag = FPU_Special(st0_ptr);
258
259 switch (tag) {
260 case TW_Denormal:
261 if (denormal_operand() < 0)
262 return;
263 goto denormal_arg;
264 case TW_Infinity:
265 if (signnegative(st0_ptr)) {
266 /* -infinity gives -1 (p16-10) */
267 FPU_copy_to_reg0(&CONST_1, TAG_Valid);
268 setnegative(st0_ptr);
269 }
270 return;
271 default:
272 single_arg_error(st0_ptr, tag);
273 }
274}
275
276static void fptan(FPU_REG *st0_ptr, u_char st0_tag)
277{
278 FPU_REG *st_new_ptr;
279 int q;
280 u_char arg_sign = getsign(st0_ptr);
281
282 /* Stack underflow has higher priority */
283 if (st0_tag == TAG_Empty) {
284 FPU_stack_underflow(); /* Puts a QNaN in st(0) */
285 if (control_word & CW_Invalid) {
286 st_new_ptr = &st(-1);
287 push();
288 FPU_stack_underflow(); /* Puts a QNaN in the new st(0) */
289 }
290 return;
291 }
292
293 if (STACK_OVERFLOW) {
294 FPU_stack_overflow();
295 return;
296 }
297
298 if (st0_tag == TAG_Valid) {
299 if (exponent(st0_ptr) > -40) {
300 if ((q = trig_arg(st0_ptr, 0)) == -1) {
301 /* Operand is out of range */
302 return;
303 }
304
305 poly_tan(st0_ptr);
306 setsign(st0_ptr, (q & 1) ^ (arg_sign != 0));
307 set_precision_flag_up(); /* We do not really know if up or down */
308 } else {
309 /* For a small arg, the result == the argument */
310 /* Underflow may happen */
311
312 denormal_arg:
313
314 FPU_to_exp16(st0_ptr, st0_ptr);
315
316 st0_tag =
317 FPU_round(st0_ptr, 1, 0, FULL_PRECISION, arg_sign);
318 FPU_settag0(st0_tag);
319 }
320 push();
321 FPU_copy_to_reg0(&CONST_1, TAG_Valid);
322 return;
323 }
324
325 if (st0_tag == TAG_Zero) {
326 push();
327 FPU_copy_to_reg0(&CONST_1, TAG_Valid);
328 setcc(0);
329 return;
330 }
331
332 if (st0_tag == TAG_Special)
333 st0_tag = FPU_Special(st0_ptr);
334
335 if (st0_tag == TW_Denormal) {
336 if (denormal_operand() < 0)
337 return;
338
339 goto denormal_arg;
340 }
341
342 if (st0_tag == TW_Infinity) {
343 /* The 80486 treats infinity as an invalid operand */
344 if (arith_invalid(0) >= 0) {
345 st_new_ptr = &st(-1);
346 push();
347 arith_invalid(0);
348 }
349 return;
350 }
351
352 single_arg_2_error(st0_ptr, st0_tag);
353}
354
355static void fxtract(FPU_REG *st0_ptr, u_char st0_tag)
356{
357 FPU_REG *st_new_ptr;
358 u_char sign;
359 register FPU_REG *st1_ptr = st0_ptr; /* anticipate */
360
361 if (STACK_OVERFLOW) {
362 FPU_stack_overflow();
363 return;
364 }
365
366 clear_C1();
367
368 if (st0_tag == TAG_Valid) {
369 long e;
370
371 push();
372 sign = getsign(st1_ptr);
373 reg_copy(st1_ptr, st_new_ptr);
374 setexponent16(st_new_ptr, exponent(st_new_ptr));
375
376 denormal_arg:
377
378 e = exponent16(st_new_ptr);
379 convert_l2reg(&e, 1);
380 setexponentpos(st_new_ptr, 0);
381 setsign(st_new_ptr, sign);
382 FPU_settag0(TAG_Valid); /* Needed if arg was a denormal */
383 return;
384 } else if (st0_tag == TAG_Zero) {
385 sign = getsign(st0_ptr);
386
387 if (FPU_divide_by_zero(0, SIGN_NEG) < 0)
388 return;
389
390 push();
391 FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
392 setsign(st_new_ptr, sign);
393 return;
394 }
395
396 if (st0_tag == TAG_Special)
397 st0_tag = FPU_Special(st0_ptr);
398
399 if (st0_tag == TW_Denormal) {
400 if (denormal_operand() < 0)
401 return;
402
403 push();
404 sign = getsign(st1_ptr);
405 FPU_to_exp16(st1_ptr, st_new_ptr);
406 goto denormal_arg;
407 } else if (st0_tag == TW_Infinity) {
408 sign = getsign(st0_ptr);
409 setpositive(st0_ptr);
410 push();
411 FPU_copy_to_reg0(&CONST_INF, TAG_Special);
412 setsign(st_new_ptr, sign);
413 return;
414 } else if (st0_tag == TW_NaN) {
415 if (real_1op_NaN(st0_ptr) < 0)
416 return;
417
418 push();
419 FPU_copy_to_reg0(st0_ptr, TAG_Special);
420 return;
421 } else if (st0_tag == TAG_Empty) {
422 /* Is this the correct behaviour? */
423 if (control_word & EX_Invalid) {
424 FPU_stack_underflow();
425 push();
426 FPU_stack_underflow();
427 } else
428 EXCEPTION(EX_StackUnder);
429 }
430#ifdef PARANOID
431 else
432 EXCEPTION(EX_INTERNAL | 0x119);
433#endif /* PARANOID */
434}
435
436static void fdecstp(void)
437{
438 clear_C1();
439 top--;
440}
441
442static void fincstp(void)
443{
444 clear_C1();
445 top++;
446}
447
448static void fsqrt_(FPU_REG *st0_ptr, u_char st0_tag)
449{
450 int expon;
451
452 clear_C1();
453
454 if (st0_tag == TAG_Valid) {
455 u_char tag;
456
457 if (signnegative(st0_ptr)) {
458 arith_invalid(0); /* sqrt(negative) is invalid */
459 return;
460 }
461
462 /* make st(0) in [1.0 .. 4.0) */
463 expon = exponent(st0_ptr);
464
465 denormal_arg:
466
467 setexponent16(st0_ptr, (expon & 1));
468
469 /* Do the computation, the sign of the result will be positive. */
470 tag = wm_sqrt(st0_ptr, 0, 0, control_word, SIGN_POS);
471 addexponent(st0_ptr, expon >> 1);
472 FPU_settag0(tag);
473 return;
474 }
475
476 if (st0_tag == TAG_Zero)
477 return;
478
479 if (st0_tag == TAG_Special)
480 st0_tag = FPU_Special(st0_ptr);
481
482 if (st0_tag == TW_Infinity) {
483 if (signnegative(st0_ptr))
484 arith_invalid(0); /* sqrt(-Infinity) is invalid */
485 return;
486 } else if (st0_tag == TW_Denormal) {
487 if (signnegative(st0_ptr)) {
488 arith_invalid(0); /* sqrt(negative) is invalid */
489 return;
490 }
491
492 if (denormal_operand() < 0)
493 return;
494
495 FPU_to_exp16(st0_ptr, st0_ptr);
496
497 expon = exponent16(st0_ptr);
498
499 goto denormal_arg;
500 }
501
502 single_arg_error(st0_ptr, st0_tag);
503
504}
505
506static void frndint_(FPU_REG *st0_ptr, u_char st0_tag)
507{
508 int flags, tag;
509
510 if (st0_tag == TAG_Valid) {
511 u_char sign;
512
513 denormal_arg:
514
515 sign = getsign(st0_ptr);
516
517 if (exponent(st0_ptr) > 63)
518 return;
519
520 if (st0_tag == TW_Denormal) {
521 if (denormal_operand() < 0)
522 return;
523 }
524
525 /* Fortunately, this can't overflow to 2^64 */
526 if ((flags = FPU_round_to_int(st0_ptr, st0_tag)))
527 set_precision_flag(flags);
528
529 setexponent16(st0_ptr, 63);
530 tag = FPU_normalize(st0_ptr);
531 setsign(st0_ptr, sign);
532 FPU_settag0(tag);
533 return;
534 }
535
536 if (st0_tag == TAG_Zero)
537 return;
538
539 if (st0_tag == TAG_Special)
540 st0_tag = FPU_Special(st0_ptr);
541
542 if (st0_tag == TW_Denormal)
543 goto denormal_arg;
544 else if (st0_tag == TW_Infinity)
545 return;
546 else
547 single_arg_error(st0_ptr, st0_tag);
548}
549
550static int f_sin(FPU_REG *st0_ptr, u_char tag)
551{
552 u_char arg_sign = getsign(st0_ptr);
553
554 if (tag == TAG_Valid) {
555 int q;
556
557 if (exponent(st0_ptr) > -40) {
558 if ((q = trig_arg(st0_ptr, 0)) == -1) {
559 /* Operand is out of range */
560 return 1;
561 }
562
563 poly_sine(st0_ptr);
564
565 if (q & 2)
566 changesign(st0_ptr);
567
568 setsign(st0_ptr, getsign(st0_ptr) ^ arg_sign);
569
570 /* We do not really know if up or down */
571 set_precision_flag_up();
572 return 0;
573 } else {
574 /* For a small arg, the result == the argument */
575 set_precision_flag_up(); /* Must be up. */
576 return 0;
577 }
578 }
579
580 if (tag == TAG_Zero) {
581 setcc(0);
582 return 0;
583 }
584
585 if (tag == TAG_Special)
586 tag = FPU_Special(st0_ptr);
587
588 if (tag == TW_Denormal) {
589 if (denormal_operand() < 0)
590 return 1;
591
592 /* For a small arg, the result == the argument */
593 /* Underflow may happen */
594 FPU_to_exp16(st0_ptr, st0_ptr);
595
596 tag = FPU_round(st0_ptr, 1, 0, FULL_PRECISION, arg_sign);
597
598 FPU_settag0(tag);
599
600 return 0;
601 } else if (tag == TW_Infinity) {
602 /* The 80486 treats infinity as an invalid operand */
603 arith_invalid(0);
604 return 1;
605 } else {
606 single_arg_error(st0_ptr, tag);
607 return 1;
608 }
609}
610
611static void fsin(FPU_REG *st0_ptr, u_char tag)
612{
613 f_sin(st0_ptr, tag);
614}
615
616static int f_cos(FPU_REG *st0_ptr, u_char tag)
617{
618 u_char st0_sign;
619
620 st0_sign = getsign(st0_ptr);
621
622 if (tag == TAG_Valid) {
623 int q;
624
625 if (exponent(st0_ptr) > -40) {
626 if ((exponent(st0_ptr) < 0)
627 || ((exponent(st0_ptr) == 0)
628 && (significand(st0_ptr) <=
629 0xc90fdaa22168c234LL))) {
630 poly_cos(st0_ptr);
631
632 /* We do not really know if up or down */
633 set_precision_flag_down();
634
635 return 0;
636 } else if ((q = trig_arg(st0_ptr, FCOS)) != -1) {
637 poly_sine(st0_ptr);
638
639 if ((q + 1) & 2)
640 changesign(st0_ptr);
641
642 /* We do not really know if up or down */
643 set_precision_flag_down();
644
645 return 0;
646 } else {
647 /* Operand is out of range */
648 return 1;
649 }
650 } else {
651 denormal_arg:
652
653 setcc(0);
654 FPU_copy_to_reg0(&CONST_1, TAG_Valid);
655#ifdef PECULIAR_486
656 set_precision_flag_down(); /* 80486 appears to do this. */
657#else
658 set_precision_flag_up(); /* Must be up. */
659#endif /* PECULIAR_486 */
660 return 0;
661 }
662 } else if (tag == TAG_Zero) {
663 FPU_copy_to_reg0(&CONST_1, TAG_Valid);
664 setcc(0);
665 return 0;
666 }
667
668 if (tag == TAG_Special)
669 tag = FPU_Special(st0_ptr);
670
671 if (tag == TW_Denormal) {
672 if (denormal_operand() < 0)
673 return 1;
674
675 goto denormal_arg;
676 } else if (tag == TW_Infinity) {
677 /* The 80486 treats infinity as an invalid operand */
678 arith_invalid(0);
679 return 1;
680 } else {
681 single_arg_error(st0_ptr, tag); /* requires st0_ptr == &st(0) */
682 return 1;
683 }
684}
685
686static void fcos(FPU_REG *st0_ptr, u_char st0_tag)
687{
688 f_cos(st0_ptr, st0_tag);
689}
690
691static void fsincos(FPU_REG *st0_ptr, u_char st0_tag)
692{
693 FPU_REG *st_new_ptr;
694 FPU_REG arg;
695 u_char tag;
696
697 /* Stack underflow has higher priority */
698 if (st0_tag == TAG_Empty) {
699 FPU_stack_underflow(); /* Puts a QNaN in st(0) */
700 if (control_word & CW_Invalid) {
701 st_new_ptr = &st(-1);
702 push();
703 FPU_stack_underflow(); /* Puts a QNaN in the new st(0) */
704 }
705 return;
706 }
707
708 if (STACK_OVERFLOW) {
709 FPU_stack_overflow();
710 return;
711 }
712
713 if (st0_tag == TAG_Special)
714 tag = FPU_Special(st0_ptr);
715 else
716 tag = st0_tag;
717
718 if (tag == TW_NaN) {
719 single_arg_2_error(st0_ptr, TW_NaN);
720 return;
721 } else if (tag == TW_Infinity) {
722 /* The 80486 treats infinity as an invalid operand */
723 if (arith_invalid(0) >= 0) {
724 /* Masked response */
725 push();
726 arith_invalid(0);
727 }
728 return;
729 }
730
731 reg_copy(st0_ptr, &arg);
732 if (!f_sin(st0_ptr, st0_tag)) {
733 push();
734 FPU_copy_to_reg0(&arg, st0_tag);
735 f_cos(&st(0), st0_tag);
736 } else {
737 /* An error, so restore st(0) */
738 FPU_copy_to_reg0(&arg, st0_tag);
739 }
740}
741
742/*---------------------------------------------------------------------------*/
743/* The following all require two arguments: st(0) and st(1) */
744
745/* A lean, mean kernel for the fprem instructions. This relies upon
746 the division and rounding to an integer in do_fprem giving an
747 exact result. Because of this, rem_kernel() needs to deal only with
748 the least significant 64 bits, the more significant bits of the
749 result must be zero.
750 */
751static void rem_kernel(unsigned long long st0, unsigned long long *y,
752 unsigned long long st1, unsigned long long q, int n)
753{
754 int dummy;
755 unsigned long long x;
756
757 x = st0 << n;
758
759 /* Do the required multiplication and subtraction in the one operation */
760
761 /* lsw x -= lsw st1 * lsw q */
762 asm volatile ("mull %4; subl %%eax,%0; sbbl %%edx,%1":"=m"
763 (((unsigned *)&x)[0]), "=m"(((unsigned *)&x)[1]),
764 "=a"(dummy)
765 :"2"(((unsigned *)&st1)[0]), "m"(((unsigned *)&q)[0])
766 :"%dx");
767 /* msw x -= msw st1 * lsw q */
768 asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x)[1]),
769 "=a"(dummy)
770 :"1"(((unsigned *)&st1)[1]), "m"(((unsigned *)&q)[0])
771 :"%dx");
772 /* msw x -= lsw st1 * msw q */
773 asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x)[1]),
774 "=a"(dummy)
775 :"1"(((unsigned *)&st1)[0]), "m"(((unsigned *)&q)[1])
776 :"%dx");
777
778 *y = x;
779}
780
781/* Remainder of st(0) / st(1) */
782/* This routine produces exact results, i.e. there is never any
783 rounding or truncation, etc of the result. */
784static void do_fprem(FPU_REG *st0_ptr, u_char st0_tag, int round)
785{
786 FPU_REG *st1_ptr = &st(1);
787 u_char st1_tag = FPU_gettagi(1);
788
789 if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
790 FPU_REG tmp, st0, st1;
791 u_char st0_sign, st1_sign;
792 u_char tmptag;
793 int tag;
794 int old_cw;
795 int expdif;
796 long long q;
797 unsigned short saved_status;
798 int cc;
799
800 fprem_valid:
801 /* Convert registers for internal use. */
802 st0_sign = FPU_to_exp16(st0_ptr, &st0);
803 st1_sign = FPU_to_exp16(st1_ptr, &st1);
804 expdif = exponent16(&st0) - exponent16(&st1);
805
806 old_cw = control_word;
807 cc = 0;
808
809 /* We want the status following the denorm tests, but don't want
810 the status changed by the arithmetic operations. */
811 saved_status = partial_status;
812 control_word &= ~CW_RC;
813 control_word |= RC_CHOP;
814
815 if (expdif < 64) {
816 /* This should be the most common case */
817
818 if (expdif > -2) {
819 u_char sign = st0_sign ^ st1_sign;
820 tag = FPU_u_div(&st0, &st1, &tmp,
821 PR_64_BITS | RC_CHOP | 0x3f,
822 sign);
823 setsign(&tmp, sign);
824
825 if (exponent(&tmp) >= 0) {
826 FPU_round_to_int(&tmp, tag); /* Fortunately, this can't
827 overflow to 2^64 */
828 q = significand(&tmp);
829
830 rem_kernel(significand(&st0),
831 &significand(&tmp),
832 significand(&st1),
833 q, expdif);
834
835 setexponent16(&tmp, exponent16(&st1));
836 } else {
837 reg_copy(&st0, &tmp);
838 q = 0;
839 }
840
841 if ((round == RC_RND)
842 && (tmp.sigh & 0xc0000000)) {
843 /* We may need to subtract st(1) once more,
844 to get a result <= 1/2 of st(1). */
845 unsigned long long x;
846 expdif =
847 exponent16(&st1) - exponent16(&tmp);
848 if (expdif <= 1) {
849 if (expdif == 0)
850 x = significand(&st1) -
851 significand(&tmp);
852 else /* expdif is 1 */
853 x = (significand(&st1)
854 << 1) -
855 significand(&tmp);
856 if ((x < significand(&tmp)) ||
857 /* or equi-distant (from 0 & st(1)) and q is odd */
858 ((x == significand(&tmp))
859 && (q & 1))) {
860 st0_sign = !st0_sign;
861 significand(&tmp) = x;
862 q++;
863 }
864 }
865 }
866
867 if (q & 4)
868 cc |= SW_C0;
869 if (q & 2)
870 cc |= SW_C3;
871 if (q & 1)
872 cc |= SW_C1;
873 } else {
874 control_word = old_cw;
875 setcc(0);
876 return;
877 }
878 } else {
879 /* There is a large exponent difference ( >= 64 ) */
880 /* To make much sense, the code in this section should
881 be done at high precision. */
882 int exp_1, N;
883 u_char sign;
884
885 /* prevent overflow here */
886 /* N is 'a number between 32 and 63' (p26-113) */
887 reg_copy(&st0, &tmp);
888 tmptag = st0_tag;
889 N = (expdif & 0x0000001f) + 32; /* This choice gives results
890 identical to an AMD 486 */
891 setexponent16(&tmp, N);
892 exp_1 = exponent16(&st1);
893 setexponent16(&st1, 0);
894 expdif -= N;
895
896 sign = getsign(&tmp) ^ st1_sign;
897 tag =
898 FPU_u_div(&tmp, &st1, &tmp,
899 PR_64_BITS | RC_CHOP | 0x3f, sign);
900 setsign(&tmp, sign);
901
902 FPU_round_to_int(&tmp, tag); /* Fortunately, this can't
903 overflow to 2^64 */
904
905 rem_kernel(significand(&st0),
906 &significand(&tmp),
907 significand(&st1),
908 significand(&tmp), exponent(&tmp)
909 );
910 setexponent16(&tmp, exp_1 + expdif);
911
912 /* It is possible for the operation to be complete here.
913 What does the IEEE standard say? The Intel 80486 manual
914 implies that the operation will never be completed at this
915 point, and the behaviour of a real 80486 confirms this.
916 */
917 if (!(tmp.sigh | tmp.sigl)) {
918 /* The result is zero */
919 control_word = old_cw;
920 partial_status = saved_status;
921 FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
922 setsign(&st0, st0_sign);
923#ifdef PECULIAR_486
924 setcc(SW_C2);
925#else
926 setcc(0);
927#endif /* PECULIAR_486 */
928 return;
929 }
930 cc = SW_C2;
931 }
932
933 control_word = old_cw;
934 partial_status = saved_status;
935 tag = FPU_normalize_nuo(&tmp);
936 reg_copy(&tmp, st0_ptr);
937
938 /* The only condition to be looked for is underflow,
939 and it can occur here only if underflow is unmasked. */
940 if ((exponent16(&tmp) <= EXP_UNDER) && (tag != TAG_Zero)
941 && !(control_word & CW_Underflow)) {
942 setcc(cc);
943 tag = arith_underflow(st0_ptr);
944 setsign(st0_ptr, st0_sign);
945 FPU_settag0(tag);
946 return;
947 } else if ((exponent16(&tmp) > EXP_UNDER) || (tag == TAG_Zero)) {
948 stdexp(st0_ptr);
949 setsign(st0_ptr, st0_sign);
950 } else {
951 tag =
952 FPU_round(st0_ptr, 0, 0, FULL_PRECISION, st0_sign);
953 }
954 FPU_settag0(tag);
955 setcc(cc);
956
957 return;
958 }
959
960 if (st0_tag == TAG_Special)
961 st0_tag = FPU_Special(st0_ptr);
962 if (st1_tag == TAG_Special)
963 st1_tag = FPU_Special(st1_ptr);
964
965 if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
966 || ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
967 || ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
968 if (denormal_operand() < 0)
969 return;
970 goto fprem_valid;
971 } else if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
972 FPU_stack_underflow();
973 return;
974 } else if (st0_tag == TAG_Zero) {
975 if (st1_tag == TAG_Valid) {
976 setcc(0);
977 return;
978 } else if (st1_tag == TW_Denormal) {
979 if (denormal_operand() < 0)
980 return;
981 setcc(0);
982 return;
983 } else if (st1_tag == TAG_Zero) {
984 arith_invalid(0);
985 return;
986 } /* fprem(?,0) always invalid */
987 else if (st1_tag == TW_Infinity) {
988 setcc(0);
989 return;
990 }
991 } else if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
992 if (st1_tag == TAG_Zero) {
993 arith_invalid(0); /* fprem(Valid,Zero) is invalid */
994 return;
995 } else if (st1_tag != TW_NaN) {
996 if (((st0_tag == TW_Denormal)
997 || (st1_tag == TW_Denormal))
998 && (denormal_operand() < 0))
999 return;
1000
1001 if (st1_tag == TW_Infinity) {
1002 /* fprem(Valid,Infinity) is o.k. */
1003 setcc(0);
1004 return;
1005 }
1006 }
1007 } else if (st0_tag == TW_Infinity) {
1008 if (st1_tag != TW_NaN) {
1009 arith_invalid(0); /* fprem(Infinity,?) is invalid */
1010 return;
1011 }
1012 }
1013
1014 /* One of the registers must contain a NaN if we got here. */
1015
1016#ifdef PARANOID
1017 if ((st0_tag != TW_NaN) && (st1_tag != TW_NaN))
1018 EXCEPTION(EX_INTERNAL | 0x118);
1019#endif /* PARANOID */
1020
1021 real_2op_NaN(st1_ptr, st1_tag, 0, st1_ptr);
1022
1023}
1024
1025/* ST(1) <- ST(1) * log ST; pop ST */
1026static void fyl2x(FPU_REG *st0_ptr, u_char st0_tag)
1027{
1028 FPU_REG *st1_ptr = &st(1), exponent;
1029 u_char st1_tag = FPU_gettagi(1);
1030 u_char sign;
1031 int e, tag;
1032
1033 clear_C1();
1034
1035 if ((st0_tag == TAG_Valid) && (st1_tag == TAG_Valid)) {
1036 both_valid:
1037 /* Both regs are Valid or Denormal */
1038 if (signpositive(st0_ptr)) {
1039 if (st0_tag == TW_Denormal)
1040 FPU_to_exp16(st0_ptr, st0_ptr);
1041 else
1042 /* Convert st(0) for internal use. */
1043 setexponent16(st0_ptr, exponent(st0_ptr));
1044
1045 if ((st0_ptr->sigh == 0x80000000)
1046 && (st0_ptr->sigl == 0)) {
1047 /* Special case. The result can be precise. */
1048 u_char esign;
1049 e = exponent16(st0_ptr);
1050 if (e >= 0) {
1051 exponent.sigh = e;
1052 esign = SIGN_POS;
1053 } else {
1054 exponent.sigh = -e;
1055 esign = SIGN_NEG;
1056 }
1057 exponent.sigl = 0;
1058 setexponent16(&exponent, 31);
1059 tag = FPU_normalize_nuo(&exponent);
1060 stdexp(&exponent);
1061 setsign(&exponent, esign);
1062 tag =
1063 FPU_mul(&exponent, tag, 1, FULL_PRECISION);
1064 if (tag >= 0)
1065 FPU_settagi(1, tag);
1066 } else {
1067 /* The usual case */
1068 sign = getsign(st1_ptr);
1069 if (st1_tag == TW_Denormal)
1070 FPU_to_exp16(st1_ptr, st1_ptr);
1071 else
1072 /* Convert st(1) for internal use. */
1073 setexponent16(st1_ptr,
1074 exponent(st1_ptr));
1075 poly_l2(st0_ptr, st1_ptr, sign);
1076 }
1077 } else {
1078 /* negative */
1079 if (arith_invalid(1) < 0)
1080 return;
1081 }
1082
1083 FPU_pop();
1084
1085 return;
1086 }
1087
1088 if (st0_tag == TAG_Special)
1089 st0_tag = FPU_Special(st0_ptr);
1090 if (st1_tag == TAG_Special)
1091 st1_tag = FPU_Special(st1_ptr);
1092
1093 if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
1094 FPU_stack_underflow_pop(1);
1095 return;
1096 } else if ((st0_tag <= TW_Denormal) && (st1_tag <= TW_Denormal)) {
1097 if (st0_tag == TAG_Zero) {
1098 if (st1_tag == TAG_Zero) {
1099 /* Both args zero is invalid */
1100 if (arith_invalid(1) < 0)
1101 return;
1102 } else {
1103 u_char sign;
1104 sign = getsign(st1_ptr) ^ SIGN_NEG;
1105 if (FPU_divide_by_zero(1, sign) < 0)
1106 return;
1107
1108 setsign(st1_ptr, sign);
1109 }
1110 } else if (st1_tag == TAG_Zero) {
1111 /* st(1) contains zero, st(0) valid <> 0 */
1112 /* Zero is the valid answer */
1113 sign = getsign(st1_ptr);
1114
1115 if (signnegative(st0_ptr)) {
1116 /* log(negative) */
1117 if (arith_invalid(1) < 0)
1118 return;
1119 } else if ((st0_tag == TW_Denormal)
1120 && (denormal_operand() < 0))
1121 return;
1122 else {
1123 if (exponent(st0_ptr) < 0)
1124 sign ^= SIGN_NEG;
1125
1126 FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
1127 setsign(st1_ptr, sign);
1128 }
1129 } else {
1130 /* One or both operands are denormals. */
1131 if (denormal_operand() < 0)
1132 return;
1133 goto both_valid;
1134 }
1135 } else if ((st0_tag == TW_NaN) || (st1_tag == TW_NaN)) {
1136 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1137 return;
1138 }
1139 /* One or both arg must be an infinity */
1140 else if (st0_tag == TW_Infinity) {
1141 if ((signnegative(st0_ptr)) || (st1_tag == TAG_Zero)) {
1142 /* log(-infinity) or 0*log(infinity) */
1143 if (arith_invalid(1) < 0)
1144 return;
1145 } else {
1146 u_char sign = getsign(st1_ptr);
1147
1148 if ((st1_tag == TW_Denormal)
1149 && (denormal_operand() < 0))
1150 return;
1151
1152 FPU_copy_to_reg1(&CONST_INF, TAG_Special);
1153 setsign(st1_ptr, sign);
1154 }
1155 }
1156 /* st(1) must be infinity here */
1157 else if (((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal))
1158 && (signpositive(st0_ptr))) {
1159 if (exponent(st0_ptr) >= 0) {
1160 if ((exponent(st0_ptr) == 0) &&
1161 (st0_ptr->sigh == 0x80000000) &&
1162 (st0_ptr->sigl == 0)) {
1163 /* st(0) holds 1.0 */
1164 /* infinity*log(1) */
1165 if (arith_invalid(1) < 0)
1166 return;
1167 }
1168 /* else st(0) is positive and > 1.0 */
1169 } else {
1170 /* st(0) is positive and < 1.0 */
1171
1172 if ((st0_tag == TW_Denormal)
1173 && (denormal_operand() < 0))
1174 return;
1175
1176 changesign(st1_ptr);
1177 }
1178 } else {
1179 /* st(0) must be zero or negative */
1180 if (st0_tag == TAG_Zero) {
1181 /* This should be invalid, but a real 80486 is happy with it. */
1182
1183#ifndef PECULIAR_486
1184 sign = getsign(st1_ptr);
1185 if (FPU_divide_by_zero(1, sign) < 0)
1186 return;
1187#endif /* PECULIAR_486 */
1188
1189 changesign(st1_ptr);
1190 } else if (arith_invalid(1) < 0) /* log(negative) */
1191 return;
1192 }
1193
1194 FPU_pop();
1195}
1196
1197static void fpatan(FPU_REG *st0_ptr, u_char st0_tag)
1198{
1199 FPU_REG *st1_ptr = &st(1);
1200 u_char st1_tag = FPU_gettagi(1);
1201 int tag;
1202
1203 clear_C1();
1204 if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1205 valid_atan:
1206
1207 poly_atan(st0_ptr, st0_tag, st1_ptr, st1_tag);
1208
1209 FPU_pop();
1210
1211 return;
1212 }
1213
1214 if (st0_tag == TAG_Special)
1215 st0_tag = FPU_Special(st0_ptr);
1216 if (st1_tag == TAG_Special)
1217 st1_tag = FPU_Special(st1_ptr);
1218
1219 if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
1220 || ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
1221 || ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
1222 if (denormal_operand() < 0)
1223 return;
1224
1225 goto valid_atan;
1226 } else if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
1227 FPU_stack_underflow_pop(1);
1228 return;
1229 } else if ((st0_tag == TW_NaN) || (st1_tag == TW_NaN)) {
1230 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) >= 0)
1231 FPU_pop();
1232 return;
1233 } else if ((st0_tag == TW_Infinity) || (st1_tag == TW_Infinity)) {
1234 u_char sign = getsign(st1_ptr);
1235 if (st0_tag == TW_Infinity) {
1236 if (st1_tag == TW_Infinity) {
1237 if (signpositive(st0_ptr)) {
1238 FPU_copy_to_reg1(&CONST_PI4, TAG_Valid);
1239 } else {
1240 setpositive(st1_ptr);
1241 tag =
1242 FPU_u_add(&CONST_PI4, &CONST_PI2,
1243 st1_ptr, FULL_PRECISION,
1244 SIGN_POS,
1245 exponent(&CONST_PI4),
1246 exponent(&CONST_PI2));
1247 if (tag >= 0)
1248 FPU_settagi(1, tag);
1249 }
1250 } else {
1251 if ((st1_tag == TW_Denormal)
1252 && (denormal_operand() < 0))
1253 return;
1254
1255 if (signpositive(st0_ptr)) {
1256 FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
1257 setsign(st1_ptr, sign); /* An 80486 preserves the sign */
1258 FPU_pop();
1259 return;
1260 } else {
1261 FPU_copy_to_reg1(&CONST_PI, TAG_Valid);
1262 }
1263 }
1264 } else {
1265 /* st(1) is infinity, st(0) not infinity */
1266 if ((st0_tag == TW_Denormal)
1267 && (denormal_operand() < 0))
1268 return;
1269
1270 FPU_copy_to_reg1(&CONST_PI2, TAG_Valid);
1271 }
1272 setsign(st1_ptr, sign);
1273 } else if (st1_tag == TAG_Zero) {
1274 /* st(0) must be valid or zero */
1275 u_char sign = getsign(st1_ptr);
1276
1277 if ((st0_tag == TW_Denormal) && (denormal_operand() < 0))
1278 return;
1279
1280 if (signpositive(st0_ptr)) {
1281 /* An 80486 preserves the sign */
1282 FPU_pop();
1283 return;
1284 }
1285
1286 FPU_copy_to_reg1(&CONST_PI, TAG_Valid);
1287 setsign(st1_ptr, sign);
1288 } else if (st0_tag == TAG_Zero) {
1289 /* st(1) must be TAG_Valid here */
1290 u_char sign = getsign(st1_ptr);
1291
1292 if ((st1_tag == TW_Denormal) && (denormal_operand() < 0))
1293 return;
1294
1295 FPU_copy_to_reg1(&CONST_PI2, TAG_Valid);
1296 setsign(st1_ptr, sign);
1297 }
1298#ifdef PARANOID
1299 else
1300 EXCEPTION(EX_INTERNAL | 0x125);
1301#endif /* PARANOID */
1302
1303 FPU_pop();
1304 set_precision_flag_up(); /* We do not really know if up or down */
1305}
1306
1307static void fprem(FPU_REG *st0_ptr, u_char st0_tag)
1308{
1309 do_fprem(st0_ptr, st0_tag, RC_CHOP);
1310}
1311
1312static void fprem1(FPU_REG *st0_ptr, u_char st0_tag)
1313{
1314 do_fprem(st0_ptr, st0_tag, RC_RND);
1315}
1316
1317static void fyl2xp1(FPU_REG *st0_ptr, u_char st0_tag)
1318{
1319 u_char sign, sign1;
1320 FPU_REG *st1_ptr = &st(1), a, b;
1321 u_char st1_tag = FPU_gettagi(1);
1322
1323 clear_C1();
1324 if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1325 valid_yl2xp1:
1326
1327 sign = getsign(st0_ptr);
1328 sign1 = getsign(st1_ptr);
1329
1330 FPU_to_exp16(st0_ptr, &a);
1331 FPU_to_exp16(st1_ptr, &b);
1332
1333 if (poly_l2p1(sign, sign1, &a, &b, st1_ptr))
1334 return;
1335
1336 FPU_pop();
1337 return;
1338 }
1339
1340 if (st0_tag == TAG_Special)
1341 st0_tag = FPU_Special(st0_ptr);
1342 if (st1_tag == TAG_Special)
1343 st1_tag = FPU_Special(st1_ptr);
1344
1345 if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
1346 || ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
1347 || ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
1348 if (denormal_operand() < 0)
1349 return;
1350
1351 goto valid_yl2xp1;
1352 } else if ((st0_tag == TAG_Empty) | (st1_tag == TAG_Empty)) {
1353 FPU_stack_underflow_pop(1);
1354 return;
1355 } else if (st0_tag == TAG_Zero) {
1356 switch (st1_tag) {
1357 case TW_Denormal:
1358 if (denormal_operand() < 0)
1359 return;
1360 fallthrough;
1361 case TAG_Zero:
1362 case TAG_Valid:
1363 setsign(st0_ptr, getsign(st0_ptr) ^ getsign(st1_ptr));
1364 FPU_copy_to_reg1(st0_ptr, st0_tag);
1365 break;
1366
1367 case TW_Infinity:
1368 /* Infinity*log(1) */
1369 if (arith_invalid(1) < 0)
1370 return;
1371 break;
1372
1373 case TW_NaN:
1374 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1375 return;
1376 break;
1377
1378 default:
1379#ifdef PARANOID
1380 EXCEPTION(EX_INTERNAL | 0x116);
1381 return;
1382#endif /* PARANOID */
1383 break;
1384 }
1385 } else if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
1386 switch (st1_tag) {
1387 case TAG_Zero:
1388 if (signnegative(st0_ptr)) {
1389 if (exponent(st0_ptr) >= 0) {
1390 /* st(0) holds <= -1.0 */
1391#ifdef PECULIAR_486 /* Stupid 80486 doesn't worry about log(negative). */
1392 changesign(st1_ptr);
1393#else
1394 if (arith_invalid(1) < 0)
1395 return;
1396#endif /* PECULIAR_486 */
1397 } else if ((st0_tag == TW_Denormal)
1398 && (denormal_operand() < 0))
1399 return;
1400 else
1401 changesign(st1_ptr);
1402 } else if ((st0_tag == TW_Denormal)
1403 && (denormal_operand() < 0))
1404 return;
1405 break;
1406
1407 case TW_Infinity:
1408 if (signnegative(st0_ptr)) {
1409 if ((exponent(st0_ptr) >= 0) &&
1410 !((st0_ptr->sigh == 0x80000000) &&
1411 (st0_ptr->sigl == 0))) {
1412 /* st(0) holds < -1.0 */
1413#ifdef PECULIAR_486 /* Stupid 80486 doesn't worry about log(negative). */
1414 changesign(st1_ptr);
1415#else
1416 if (arith_invalid(1) < 0)
1417 return;
1418#endif /* PECULIAR_486 */
1419 } else if ((st0_tag == TW_Denormal)
1420 && (denormal_operand() < 0))
1421 return;
1422 else
1423 changesign(st1_ptr);
1424 } else if ((st0_tag == TW_Denormal)
1425 && (denormal_operand() < 0))
1426 return;
1427 break;
1428
1429 case TW_NaN:
1430 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1431 return;
1432 }
1433
1434 } else if (st0_tag == TW_NaN) {
1435 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1436 return;
1437 } else if (st0_tag == TW_Infinity) {
1438 if (st1_tag == TW_NaN) {
1439 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1440 return;
1441 } else if (signnegative(st0_ptr)) {
1442#ifndef PECULIAR_486
1443 /* This should have higher priority than denormals, but... */
1444 if (arith_invalid(1) < 0) /* log(-infinity) */
1445 return;
1446#endif /* PECULIAR_486 */
1447 if ((st1_tag == TW_Denormal)
1448 && (denormal_operand() < 0))
1449 return;
1450#ifdef PECULIAR_486
1451 /* Denormal operands actually get higher priority */
1452 if (arith_invalid(1) < 0) /* log(-infinity) */
1453 return;
1454#endif /* PECULIAR_486 */
1455 } else if (st1_tag == TAG_Zero) {
1456 /* log(infinity) */
1457 if (arith_invalid(1) < 0)
1458 return;
1459 }
1460
1461 /* st(1) must be valid here. */
1462
1463 else if ((st1_tag == TW_Denormal) && (denormal_operand() < 0))
1464 return;
1465
1466 /* The Manual says that log(Infinity) is invalid, but a real
1467 80486 sensibly says that it is o.k. */
1468 else {
1469 u_char sign = getsign(st1_ptr);
1470 FPU_copy_to_reg1(&CONST_INF, TAG_Special);
1471 setsign(st1_ptr, sign);
1472 }
1473 }
1474#ifdef PARANOID
1475 else {
1476 EXCEPTION(EX_INTERNAL | 0x117);
1477 return;
1478 }
1479#endif /* PARANOID */
1480
1481 FPU_pop();
1482 return;
1483
1484}
1485
1486static void fscale(FPU_REG *st0_ptr, u_char st0_tag)
1487{
1488 FPU_REG *st1_ptr = &st(1);
1489 u_char st1_tag = FPU_gettagi(1);
1490 int old_cw = control_word;
1491 u_char sign = getsign(st0_ptr);
1492
1493 clear_C1();
1494 if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1495 long scale;
1496 FPU_REG tmp;
1497
1498 /* Convert register for internal use. */
1499 setexponent16(st0_ptr, exponent(st0_ptr));
1500
1501 valid_scale:
1502
1503 if (exponent(st1_ptr) > 30) {
1504 /* 2^31 is far too large, would require 2^(2^30) or 2^(-2^30) */
1505
1506 if (signpositive(st1_ptr)) {
1507 EXCEPTION(EX_Overflow);
1508 FPU_copy_to_reg0(&CONST_INF, TAG_Special);
1509 } else {
1510 EXCEPTION(EX_Underflow);
1511 FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
1512 }
1513 setsign(st0_ptr, sign);
1514 return;
1515 }
1516
1517 control_word &= ~CW_RC;
1518 control_word |= RC_CHOP;
1519 reg_copy(st1_ptr, &tmp);
1520 FPU_round_to_int(&tmp, st1_tag); /* This can never overflow here */
1521 control_word = old_cw;
1522 scale = signnegative(st1_ptr) ? -tmp.sigl : tmp.sigl;
1523 scale += exponent16(st0_ptr);
1524
1525 setexponent16(st0_ptr, scale);
1526
1527 /* Use FPU_round() to properly detect under/overflow etc */
1528 FPU_round(st0_ptr, 0, 0, control_word, sign);
1529
1530 return;
1531 }
1532
1533 if (st0_tag == TAG_Special)
1534 st0_tag = FPU_Special(st0_ptr);
1535 if (st1_tag == TAG_Special)
1536 st1_tag = FPU_Special(st1_ptr);
1537
1538 if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
1539 switch (st1_tag) {
1540 case TAG_Valid:
1541 /* st(0) must be a denormal */
1542 if ((st0_tag == TW_Denormal)
1543 && (denormal_operand() < 0))
1544 return;
1545
1546 FPU_to_exp16(st0_ptr, st0_ptr); /* Will not be left on stack */
1547 goto valid_scale;
1548
1549 case TAG_Zero:
1550 if (st0_tag == TW_Denormal)
1551 denormal_operand();
1552 return;
1553
1554 case TW_Denormal:
1555 denormal_operand();
1556 return;
1557
1558 case TW_Infinity:
1559 if ((st0_tag == TW_Denormal)
1560 && (denormal_operand() < 0))
1561 return;
1562
1563 if (signpositive(st1_ptr))
1564 FPU_copy_to_reg0(&CONST_INF, TAG_Special);
1565 else
1566 FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
1567 setsign(st0_ptr, sign);
1568 return;
1569
1570 case TW_NaN:
1571 real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1572 return;
1573 }
1574 } else if (st0_tag == TAG_Zero) {
1575 switch (st1_tag) {
1576 case TAG_Valid:
1577 case TAG_Zero:
1578 return;
1579
1580 case TW_Denormal:
1581 denormal_operand();
1582 return;
1583
1584 case TW_Infinity:
1585 if (signpositive(st1_ptr))
1586 arith_invalid(0); /* Zero scaled by +Infinity */
1587 return;
1588
1589 case TW_NaN:
1590 real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1591 return;
1592 }
1593 } else if (st0_tag == TW_Infinity) {
1594 switch (st1_tag) {
1595 case TAG_Valid:
1596 case TAG_Zero:
1597 return;
1598
1599 case TW_Denormal:
1600 denormal_operand();
1601 return;
1602
1603 case TW_Infinity:
1604 if (signnegative(st1_ptr))
1605 arith_invalid(0); /* Infinity scaled by -Infinity */
1606 return;
1607
1608 case TW_NaN:
1609 real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1610 return;
1611 }
1612 } else if (st0_tag == TW_NaN) {
1613 if (st1_tag != TAG_Empty) {
1614 real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1615 return;
1616 }
1617 }
1618#ifdef PARANOID
1619 if (!((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty))) {
1620 EXCEPTION(EX_INTERNAL | 0x115);
1621 return;
1622 }
1623#endif
1624
1625 /* At least one of st(0), st(1) must be empty */
1626 FPU_stack_underflow();
1627
1628}
1629
1630/*---------------------------------------------------------------------------*/
1631
1632static FUNC_ST0 const trig_table_a[] = {
1633 f2xm1, fyl2x, fptan, fpatan,
1634 fxtract, fprem1, (FUNC_ST0) fdecstp, (FUNC_ST0) fincstp
1635};
1636
1637void FPU_triga(void)
1638{
1639 (trig_table_a[FPU_rm]) (&st(0), FPU_gettag0());
1640}
1641
1642static FUNC_ST0 const trig_table_b[] = {
1643 fprem, fyl2xp1, fsqrt_, fsincos, frndint_, fscale, fsin, fcos
1644};
1645
1646void FPU_trigb(void)
1647{
1648 (trig_table_b[FPU_rm]) (&st(0), FPU_gettag0());
1649}
1/*---------------------------------------------------------------------------+
2 | fpu_trig.c |
3 | |
4 | Implementation of the FPU "transcendental" functions. |
5 | |
6 | Copyright (C) 1992,1993,1994,1997,1999 |
7 | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, |
8 | Australia. E-mail billm@melbpc.org.au |
9 | |
10 | |
11 +---------------------------------------------------------------------------*/
12
13#include "fpu_system.h"
14#include "exception.h"
15#include "fpu_emu.h"
16#include "status_w.h"
17#include "control_w.h"
18#include "reg_constant.h"
19
20static void rem_kernel(unsigned long long st0, unsigned long long *y,
21 unsigned long long st1, unsigned long long q, int n);
22
23#define BETTER_THAN_486
24
25#define FCOS 4
26
27/* Used only by fptan, fsin, fcos, and fsincos. */
28/* This routine produces very accurate results, similar to
29 using a value of pi with more than 128 bits precision. */
30/* Limited measurements show no results worse than 64 bit precision
31 except for the results for arguments close to 2^63, where the
32 precision of the result sometimes degrades to about 63.9 bits */
33static int trig_arg(FPU_REG *st0_ptr, int even)
34{
35 FPU_REG tmp;
36 u_char tmptag;
37 unsigned long long q;
38 int old_cw = control_word, saved_status = partial_status;
39 int tag, st0_tag = TAG_Valid;
40
41 if (exponent(st0_ptr) >= 63) {
42 partial_status |= SW_C2; /* Reduction incomplete. */
43 return -1;
44 }
45
46 control_word &= ~CW_RC;
47 control_word |= RC_CHOP;
48
49 setpositive(st0_ptr);
50 tag = FPU_u_div(st0_ptr, &CONST_PI2, &tmp, PR_64_BITS | RC_CHOP | 0x3f,
51 SIGN_POS);
52
53 FPU_round_to_int(&tmp, tag); /* Fortunately, this can't overflow
54 to 2^64 */
55 q = significand(&tmp);
56 if (q) {
57 rem_kernel(significand(st0_ptr),
58 &significand(&tmp),
59 significand(&CONST_PI2),
60 q, exponent(st0_ptr) - exponent(&CONST_PI2));
61 setexponent16(&tmp, exponent(&CONST_PI2));
62 st0_tag = FPU_normalize(&tmp);
63 FPU_copy_to_reg0(&tmp, st0_tag);
64 }
65
66 if ((even && !(q & 1)) || (!even && (q & 1))) {
67 st0_tag =
68 FPU_sub(REV | LOADED | TAG_Valid, (int)&CONST_PI2,
69 FULL_PRECISION);
70
71#ifdef BETTER_THAN_486
72 /* So far, the results are exact but based upon a 64 bit
73 precision approximation to pi/2. The technique used
74 now is equivalent to using an approximation to pi/2 which
75 is accurate to about 128 bits. */
76 if ((exponent(st0_ptr) <= exponent(&CONST_PI2extra) + 64)
77 || (q > 1)) {
78 /* This code gives the effect of having pi/2 to better than
79 128 bits precision. */
80
81 significand(&tmp) = q + 1;
82 setexponent16(&tmp, 63);
83 FPU_normalize(&tmp);
84 tmptag =
85 FPU_u_mul(&CONST_PI2extra, &tmp, &tmp,
86 FULL_PRECISION, SIGN_POS,
87 exponent(&CONST_PI2extra) +
88 exponent(&tmp));
89 setsign(&tmp, getsign(&CONST_PI2extra));
90 st0_tag = FPU_add(&tmp, tmptag, 0, FULL_PRECISION);
91 if (signnegative(st0_ptr)) {
92 /* CONST_PI2extra is negative, so the result of the addition
93 can be negative. This means that the argument is actually
94 in a different quadrant. The correction is always < pi/2,
95 so it can't overflow into yet another quadrant. */
96 setpositive(st0_ptr);
97 q++;
98 }
99 }
100#endif /* BETTER_THAN_486 */
101 }
102#ifdef BETTER_THAN_486
103 else {
104 /* So far, the results are exact but based upon a 64 bit
105 precision approximation to pi/2. The technique used
106 now is equivalent to using an approximation to pi/2 which
107 is accurate to about 128 bits. */
108 if (((q > 0)
109 && (exponent(st0_ptr) <= exponent(&CONST_PI2extra) + 64))
110 || (q > 1)) {
111 /* This code gives the effect of having p/2 to better than
112 128 bits precision. */
113
114 significand(&tmp) = q;
115 setexponent16(&tmp, 63);
116 FPU_normalize(&tmp); /* This must return TAG_Valid */
117 tmptag =
118 FPU_u_mul(&CONST_PI2extra, &tmp, &tmp,
119 FULL_PRECISION, SIGN_POS,
120 exponent(&CONST_PI2extra) +
121 exponent(&tmp));
122 setsign(&tmp, getsign(&CONST_PI2extra));
123 st0_tag = FPU_sub(LOADED | (tmptag & 0x0f), (int)&tmp,
124 FULL_PRECISION);
125 if ((exponent(st0_ptr) == exponent(&CONST_PI2)) &&
126 ((st0_ptr->sigh > CONST_PI2.sigh)
127 || ((st0_ptr->sigh == CONST_PI2.sigh)
128 && (st0_ptr->sigl > CONST_PI2.sigl)))) {
129 /* CONST_PI2extra is negative, so the result of the
130 subtraction can be larger than pi/2. This means
131 that the argument is actually in a different quadrant.
132 The correction is always < pi/2, so it can't overflow
133 into yet another quadrant. */
134 st0_tag =
135 FPU_sub(REV | LOADED | TAG_Valid,
136 (int)&CONST_PI2, FULL_PRECISION);
137 q++;
138 }
139 }
140 }
141#endif /* BETTER_THAN_486 */
142
143 FPU_settag0(st0_tag);
144 control_word = old_cw;
145 partial_status = saved_status & ~SW_C2; /* Reduction complete. */
146
147 return (q & 3) | even;
148}
149
150/* Convert a long to register */
151static void convert_l2reg(long const *arg, int deststnr)
152{
153 int tag;
154 long num = *arg;
155 u_char sign;
156 FPU_REG *dest = &st(deststnr);
157
158 if (num == 0) {
159 FPU_copy_to_regi(&CONST_Z, TAG_Zero, deststnr);
160 return;
161 }
162
163 if (num > 0) {
164 sign = SIGN_POS;
165 } else {
166 num = -num;
167 sign = SIGN_NEG;
168 }
169
170 dest->sigh = num;
171 dest->sigl = 0;
172 setexponent16(dest, 31);
173 tag = FPU_normalize(dest);
174 FPU_settagi(deststnr, tag);
175 setsign(dest, sign);
176 return;
177}
178
179static void single_arg_error(FPU_REG *st0_ptr, u_char st0_tag)
180{
181 if (st0_tag == TAG_Empty)
182 FPU_stack_underflow(); /* Puts a QNaN in st(0) */
183 else if (st0_tag == TW_NaN)
184 real_1op_NaN(st0_ptr); /* return with a NaN in st(0) */
185#ifdef PARANOID
186 else
187 EXCEPTION(EX_INTERNAL | 0x0112);
188#endif /* PARANOID */
189}
190
191static void single_arg_2_error(FPU_REG *st0_ptr, u_char st0_tag)
192{
193 int isNaN;
194
195 switch (st0_tag) {
196 case TW_NaN:
197 isNaN = (exponent(st0_ptr) == EXP_OVER)
198 && (st0_ptr->sigh & 0x80000000);
199 if (isNaN && !(st0_ptr->sigh & 0x40000000)) { /* Signaling ? */
200 EXCEPTION(EX_Invalid);
201 if (control_word & CW_Invalid) {
202 /* The masked response */
203 /* Convert to a QNaN */
204 st0_ptr->sigh |= 0x40000000;
205 push();
206 FPU_copy_to_reg0(st0_ptr, TAG_Special);
207 }
208 } else if (isNaN) {
209 /* A QNaN */
210 push();
211 FPU_copy_to_reg0(st0_ptr, TAG_Special);
212 } else {
213 /* pseudoNaN or other unsupported */
214 EXCEPTION(EX_Invalid);
215 if (control_word & CW_Invalid) {
216 /* The masked response */
217 FPU_copy_to_reg0(&CONST_QNaN, TAG_Special);
218 push();
219 FPU_copy_to_reg0(&CONST_QNaN, TAG_Special);
220 }
221 }
222 break; /* return with a NaN in st(0) */
223#ifdef PARANOID
224 default:
225 EXCEPTION(EX_INTERNAL | 0x0112);
226#endif /* PARANOID */
227 }
228}
229
230/*---------------------------------------------------------------------------*/
231
232static void f2xm1(FPU_REG *st0_ptr, u_char tag)
233{
234 FPU_REG a;
235
236 clear_C1();
237
238 if (tag == TAG_Valid) {
239 /* For an 80486 FPU, the result is undefined if the arg is >= 1.0 */
240 if (exponent(st0_ptr) < 0) {
241 denormal_arg:
242
243 FPU_to_exp16(st0_ptr, &a);
244
245 /* poly_2xm1(x) requires 0 < st(0) < 1. */
246 poly_2xm1(getsign(st0_ptr), &a, st0_ptr);
247 }
248 set_precision_flag_up(); /* 80486 appears to always do this */
249 return;
250 }
251
252 if (tag == TAG_Zero)
253 return;
254
255 if (tag == TAG_Special)
256 tag = FPU_Special(st0_ptr);
257
258 switch (tag) {
259 case TW_Denormal:
260 if (denormal_operand() < 0)
261 return;
262 goto denormal_arg;
263 case TW_Infinity:
264 if (signnegative(st0_ptr)) {
265 /* -infinity gives -1 (p16-10) */
266 FPU_copy_to_reg0(&CONST_1, TAG_Valid);
267 setnegative(st0_ptr);
268 }
269 return;
270 default:
271 single_arg_error(st0_ptr, tag);
272 }
273}
274
275static void fptan(FPU_REG *st0_ptr, u_char st0_tag)
276{
277 FPU_REG *st_new_ptr;
278 int q;
279 u_char arg_sign = getsign(st0_ptr);
280
281 /* Stack underflow has higher priority */
282 if (st0_tag == TAG_Empty) {
283 FPU_stack_underflow(); /* Puts a QNaN in st(0) */
284 if (control_word & CW_Invalid) {
285 st_new_ptr = &st(-1);
286 push();
287 FPU_stack_underflow(); /* Puts a QNaN in the new st(0) */
288 }
289 return;
290 }
291
292 if (STACK_OVERFLOW) {
293 FPU_stack_overflow();
294 return;
295 }
296
297 if (st0_tag == TAG_Valid) {
298 if (exponent(st0_ptr) > -40) {
299 if ((q = trig_arg(st0_ptr, 0)) == -1) {
300 /* Operand is out of range */
301 return;
302 }
303
304 poly_tan(st0_ptr);
305 setsign(st0_ptr, (q & 1) ^ (arg_sign != 0));
306 set_precision_flag_up(); /* We do not really know if up or down */
307 } else {
308 /* For a small arg, the result == the argument */
309 /* Underflow may happen */
310
311 denormal_arg:
312
313 FPU_to_exp16(st0_ptr, st0_ptr);
314
315 st0_tag =
316 FPU_round(st0_ptr, 1, 0, FULL_PRECISION, arg_sign);
317 FPU_settag0(st0_tag);
318 }
319 push();
320 FPU_copy_to_reg0(&CONST_1, TAG_Valid);
321 return;
322 }
323
324 if (st0_tag == TAG_Zero) {
325 push();
326 FPU_copy_to_reg0(&CONST_1, TAG_Valid);
327 setcc(0);
328 return;
329 }
330
331 if (st0_tag == TAG_Special)
332 st0_tag = FPU_Special(st0_ptr);
333
334 if (st0_tag == TW_Denormal) {
335 if (denormal_operand() < 0)
336 return;
337
338 goto denormal_arg;
339 }
340
341 if (st0_tag == TW_Infinity) {
342 /* The 80486 treats infinity as an invalid operand */
343 if (arith_invalid(0) >= 0) {
344 st_new_ptr = &st(-1);
345 push();
346 arith_invalid(0);
347 }
348 return;
349 }
350
351 single_arg_2_error(st0_ptr, st0_tag);
352}
353
354static void fxtract(FPU_REG *st0_ptr, u_char st0_tag)
355{
356 FPU_REG *st_new_ptr;
357 u_char sign;
358 register FPU_REG *st1_ptr = st0_ptr; /* anticipate */
359
360 if (STACK_OVERFLOW) {
361 FPU_stack_overflow();
362 return;
363 }
364
365 clear_C1();
366
367 if (st0_tag == TAG_Valid) {
368 long e;
369
370 push();
371 sign = getsign(st1_ptr);
372 reg_copy(st1_ptr, st_new_ptr);
373 setexponent16(st_new_ptr, exponent(st_new_ptr));
374
375 denormal_arg:
376
377 e = exponent16(st_new_ptr);
378 convert_l2reg(&e, 1);
379 setexponentpos(st_new_ptr, 0);
380 setsign(st_new_ptr, sign);
381 FPU_settag0(TAG_Valid); /* Needed if arg was a denormal */
382 return;
383 } else if (st0_tag == TAG_Zero) {
384 sign = getsign(st0_ptr);
385
386 if (FPU_divide_by_zero(0, SIGN_NEG) < 0)
387 return;
388
389 push();
390 FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
391 setsign(st_new_ptr, sign);
392 return;
393 }
394
395 if (st0_tag == TAG_Special)
396 st0_tag = FPU_Special(st0_ptr);
397
398 if (st0_tag == TW_Denormal) {
399 if (denormal_operand() < 0)
400 return;
401
402 push();
403 sign = getsign(st1_ptr);
404 FPU_to_exp16(st1_ptr, st_new_ptr);
405 goto denormal_arg;
406 } else if (st0_tag == TW_Infinity) {
407 sign = getsign(st0_ptr);
408 setpositive(st0_ptr);
409 push();
410 FPU_copy_to_reg0(&CONST_INF, TAG_Special);
411 setsign(st_new_ptr, sign);
412 return;
413 } else if (st0_tag == TW_NaN) {
414 if (real_1op_NaN(st0_ptr) < 0)
415 return;
416
417 push();
418 FPU_copy_to_reg0(st0_ptr, TAG_Special);
419 return;
420 } else if (st0_tag == TAG_Empty) {
421 /* Is this the correct behaviour? */
422 if (control_word & EX_Invalid) {
423 FPU_stack_underflow();
424 push();
425 FPU_stack_underflow();
426 } else
427 EXCEPTION(EX_StackUnder);
428 }
429#ifdef PARANOID
430 else
431 EXCEPTION(EX_INTERNAL | 0x119);
432#endif /* PARANOID */
433}
434
435static void fdecstp(void)
436{
437 clear_C1();
438 top--;
439}
440
441static void fincstp(void)
442{
443 clear_C1();
444 top++;
445}
446
447static void fsqrt_(FPU_REG *st0_ptr, u_char st0_tag)
448{
449 int expon;
450
451 clear_C1();
452
453 if (st0_tag == TAG_Valid) {
454 u_char tag;
455
456 if (signnegative(st0_ptr)) {
457 arith_invalid(0); /* sqrt(negative) is invalid */
458 return;
459 }
460
461 /* make st(0) in [1.0 .. 4.0) */
462 expon = exponent(st0_ptr);
463
464 denormal_arg:
465
466 setexponent16(st0_ptr, (expon & 1));
467
468 /* Do the computation, the sign of the result will be positive. */
469 tag = wm_sqrt(st0_ptr, 0, 0, control_word, SIGN_POS);
470 addexponent(st0_ptr, expon >> 1);
471 FPU_settag0(tag);
472 return;
473 }
474
475 if (st0_tag == TAG_Zero)
476 return;
477
478 if (st0_tag == TAG_Special)
479 st0_tag = FPU_Special(st0_ptr);
480
481 if (st0_tag == TW_Infinity) {
482 if (signnegative(st0_ptr))
483 arith_invalid(0); /* sqrt(-Infinity) is invalid */
484 return;
485 } else if (st0_tag == TW_Denormal) {
486 if (signnegative(st0_ptr)) {
487 arith_invalid(0); /* sqrt(negative) is invalid */
488 return;
489 }
490
491 if (denormal_operand() < 0)
492 return;
493
494 FPU_to_exp16(st0_ptr, st0_ptr);
495
496 expon = exponent16(st0_ptr);
497
498 goto denormal_arg;
499 }
500
501 single_arg_error(st0_ptr, st0_tag);
502
503}
504
505static void frndint_(FPU_REG *st0_ptr, u_char st0_tag)
506{
507 int flags, tag;
508
509 if (st0_tag == TAG_Valid) {
510 u_char sign;
511
512 denormal_arg:
513
514 sign = getsign(st0_ptr);
515
516 if (exponent(st0_ptr) > 63)
517 return;
518
519 if (st0_tag == TW_Denormal) {
520 if (denormal_operand() < 0)
521 return;
522 }
523
524 /* Fortunately, this can't overflow to 2^64 */
525 if ((flags = FPU_round_to_int(st0_ptr, st0_tag)))
526 set_precision_flag(flags);
527
528 setexponent16(st0_ptr, 63);
529 tag = FPU_normalize(st0_ptr);
530 setsign(st0_ptr, sign);
531 FPU_settag0(tag);
532 return;
533 }
534
535 if (st0_tag == TAG_Zero)
536 return;
537
538 if (st0_tag == TAG_Special)
539 st0_tag = FPU_Special(st0_ptr);
540
541 if (st0_tag == TW_Denormal)
542 goto denormal_arg;
543 else if (st0_tag == TW_Infinity)
544 return;
545 else
546 single_arg_error(st0_ptr, st0_tag);
547}
548
549static int fsin(FPU_REG *st0_ptr, u_char tag)
550{
551 u_char arg_sign = getsign(st0_ptr);
552
553 if (tag == TAG_Valid) {
554 int q;
555
556 if (exponent(st0_ptr) > -40) {
557 if ((q = trig_arg(st0_ptr, 0)) == -1) {
558 /* Operand is out of range */
559 return 1;
560 }
561
562 poly_sine(st0_ptr);
563
564 if (q & 2)
565 changesign(st0_ptr);
566
567 setsign(st0_ptr, getsign(st0_ptr) ^ arg_sign);
568
569 /* We do not really know if up or down */
570 set_precision_flag_up();
571 return 0;
572 } else {
573 /* For a small arg, the result == the argument */
574 set_precision_flag_up(); /* Must be up. */
575 return 0;
576 }
577 }
578
579 if (tag == TAG_Zero) {
580 setcc(0);
581 return 0;
582 }
583
584 if (tag == TAG_Special)
585 tag = FPU_Special(st0_ptr);
586
587 if (tag == TW_Denormal) {
588 if (denormal_operand() < 0)
589 return 1;
590
591 /* For a small arg, the result == the argument */
592 /* Underflow may happen */
593 FPU_to_exp16(st0_ptr, st0_ptr);
594
595 tag = FPU_round(st0_ptr, 1, 0, FULL_PRECISION, arg_sign);
596
597 FPU_settag0(tag);
598
599 return 0;
600 } else if (tag == TW_Infinity) {
601 /* The 80486 treats infinity as an invalid operand */
602 arith_invalid(0);
603 return 1;
604 } else {
605 single_arg_error(st0_ptr, tag);
606 return 1;
607 }
608}
609
610static int f_cos(FPU_REG *st0_ptr, u_char tag)
611{
612 u_char st0_sign;
613
614 st0_sign = getsign(st0_ptr);
615
616 if (tag == TAG_Valid) {
617 int q;
618
619 if (exponent(st0_ptr) > -40) {
620 if ((exponent(st0_ptr) < 0)
621 || ((exponent(st0_ptr) == 0)
622 && (significand(st0_ptr) <=
623 0xc90fdaa22168c234LL))) {
624 poly_cos(st0_ptr);
625
626 /* We do not really know if up or down */
627 set_precision_flag_down();
628
629 return 0;
630 } else if ((q = trig_arg(st0_ptr, FCOS)) != -1) {
631 poly_sine(st0_ptr);
632
633 if ((q + 1) & 2)
634 changesign(st0_ptr);
635
636 /* We do not really know if up or down */
637 set_precision_flag_down();
638
639 return 0;
640 } else {
641 /* Operand is out of range */
642 return 1;
643 }
644 } else {
645 denormal_arg:
646
647 setcc(0);
648 FPU_copy_to_reg0(&CONST_1, TAG_Valid);
649#ifdef PECULIAR_486
650 set_precision_flag_down(); /* 80486 appears to do this. */
651#else
652 set_precision_flag_up(); /* Must be up. */
653#endif /* PECULIAR_486 */
654 return 0;
655 }
656 } else if (tag == TAG_Zero) {
657 FPU_copy_to_reg0(&CONST_1, TAG_Valid);
658 setcc(0);
659 return 0;
660 }
661
662 if (tag == TAG_Special)
663 tag = FPU_Special(st0_ptr);
664
665 if (tag == TW_Denormal) {
666 if (denormal_operand() < 0)
667 return 1;
668
669 goto denormal_arg;
670 } else if (tag == TW_Infinity) {
671 /* The 80486 treats infinity as an invalid operand */
672 arith_invalid(0);
673 return 1;
674 } else {
675 single_arg_error(st0_ptr, tag); /* requires st0_ptr == &st(0) */
676 return 1;
677 }
678}
679
680static void fcos(FPU_REG *st0_ptr, u_char st0_tag)
681{
682 f_cos(st0_ptr, st0_tag);
683}
684
685static void fsincos(FPU_REG *st0_ptr, u_char st0_tag)
686{
687 FPU_REG *st_new_ptr;
688 FPU_REG arg;
689 u_char tag;
690
691 /* Stack underflow has higher priority */
692 if (st0_tag == TAG_Empty) {
693 FPU_stack_underflow(); /* Puts a QNaN in st(0) */
694 if (control_word & CW_Invalid) {
695 st_new_ptr = &st(-1);
696 push();
697 FPU_stack_underflow(); /* Puts a QNaN in the new st(0) */
698 }
699 return;
700 }
701
702 if (STACK_OVERFLOW) {
703 FPU_stack_overflow();
704 return;
705 }
706
707 if (st0_tag == TAG_Special)
708 tag = FPU_Special(st0_ptr);
709 else
710 tag = st0_tag;
711
712 if (tag == TW_NaN) {
713 single_arg_2_error(st0_ptr, TW_NaN);
714 return;
715 } else if (tag == TW_Infinity) {
716 /* The 80486 treats infinity as an invalid operand */
717 if (arith_invalid(0) >= 0) {
718 /* Masked response */
719 push();
720 arith_invalid(0);
721 }
722 return;
723 }
724
725 reg_copy(st0_ptr, &arg);
726 if (!fsin(st0_ptr, st0_tag)) {
727 push();
728 FPU_copy_to_reg0(&arg, st0_tag);
729 f_cos(&st(0), st0_tag);
730 } else {
731 /* An error, so restore st(0) */
732 FPU_copy_to_reg0(&arg, st0_tag);
733 }
734}
735
736/*---------------------------------------------------------------------------*/
737/* The following all require two arguments: st(0) and st(1) */
738
739/* A lean, mean kernel for the fprem instructions. This relies upon
740 the division and rounding to an integer in do_fprem giving an
741 exact result. Because of this, rem_kernel() needs to deal only with
742 the least significant 64 bits, the more significant bits of the
743 result must be zero.
744 */
745static void rem_kernel(unsigned long long st0, unsigned long long *y,
746 unsigned long long st1, unsigned long long q, int n)
747{
748 int dummy;
749 unsigned long long x;
750
751 x = st0 << n;
752
753 /* Do the required multiplication and subtraction in the one operation */
754
755 /* lsw x -= lsw st1 * lsw q */
756 asm volatile ("mull %4; subl %%eax,%0; sbbl %%edx,%1":"=m"
757 (((unsigned *)&x)[0]), "=m"(((unsigned *)&x)[1]),
758 "=a"(dummy)
759 :"2"(((unsigned *)&st1)[0]), "m"(((unsigned *)&q)[0])
760 :"%dx");
761 /* msw x -= msw st1 * lsw q */
762 asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x)[1]),
763 "=a"(dummy)
764 :"1"(((unsigned *)&st1)[1]), "m"(((unsigned *)&q)[0])
765 :"%dx");
766 /* msw x -= lsw st1 * msw q */
767 asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x)[1]),
768 "=a"(dummy)
769 :"1"(((unsigned *)&st1)[0]), "m"(((unsigned *)&q)[1])
770 :"%dx");
771
772 *y = x;
773}
774
775/* Remainder of st(0) / st(1) */
776/* This routine produces exact results, i.e. there is never any
777 rounding or truncation, etc of the result. */
778static void do_fprem(FPU_REG *st0_ptr, u_char st0_tag, int round)
779{
780 FPU_REG *st1_ptr = &st(1);
781 u_char st1_tag = FPU_gettagi(1);
782
783 if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
784 FPU_REG tmp, st0, st1;
785 u_char st0_sign, st1_sign;
786 u_char tmptag;
787 int tag;
788 int old_cw;
789 int expdif;
790 long long q;
791 unsigned short saved_status;
792 int cc;
793
794 fprem_valid:
795 /* Convert registers for internal use. */
796 st0_sign = FPU_to_exp16(st0_ptr, &st0);
797 st1_sign = FPU_to_exp16(st1_ptr, &st1);
798 expdif = exponent16(&st0) - exponent16(&st1);
799
800 old_cw = control_word;
801 cc = 0;
802
803 /* We want the status following the denorm tests, but don't want
804 the status changed by the arithmetic operations. */
805 saved_status = partial_status;
806 control_word &= ~CW_RC;
807 control_word |= RC_CHOP;
808
809 if (expdif < 64) {
810 /* This should be the most common case */
811
812 if (expdif > -2) {
813 u_char sign = st0_sign ^ st1_sign;
814 tag = FPU_u_div(&st0, &st1, &tmp,
815 PR_64_BITS | RC_CHOP | 0x3f,
816 sign);
817 setsign(&tmp, sign);
818
819 if (exponent(&tmp) >= 0) {
820 FPU_round_to_int(&tmp, tag); /* Fortunately, this can't
821 overflow to 2^64 */
822 q = significand(&tmp);
823
824 rem_kernel(significand(&st0),
825 &significand(&tmp),
826 significand(&st1),
827 q, expdif);
828
829 setexponent16(&tmp, exponent16(&st1));
830 } else {
831 reg_copy(&st0, &tmp);
832 q = 0;
833 }
834
835 if ((round == RC_RND)
836 && (tmp.sigh & 0xc0000000)) {
837 /* We may need to subtract st(1) once more,
838 to get a result <= 1/2 of st(1). */
839 unsigned long long x;
840 expdif =
841 exponent16(&st1) - exponent16(&tmp);
842 if (expdif <= 1) {
843 if (expdif == 0)
844 x = significand(&st1) -
845 significand(&tmp);
846 else /* expdif is 1 */
847 x = (significand(&st1)
848 << 1) -
849 significand(&tmp);
850 if ((x < significand(&tmp)) ||
851 /* or equi-distant (from 0 & st(1)) and q is odd */
852 ((x == significand(&tmp))
853 && (q & 1))) {
854 st0_sign = !st0_sign;
855 significand(&tmp) = x;
856 q++;
857 }
858 }
859 }
860
861 if (q & 4)
862 cc |= SW_C0;
863 if (q & 2)
864 cc |= SW_C3;
865 if (q & 1)
866 cc |= SW_C1;
867 } else {
868 control_word = old_cw;
869 setcc(0);
870 return;
871 }
872 } else {
873 /* There is a large exponent difference ( >= 64 ) */
874 /* To make much sense, the code in this section should
875 be done at high precision. */
876 int exp_1, N;
877 u_char sign;
878
879 /* prevent overflow here */
880 /* N is 'a number between 32 and 63' (p26-113) */
881 reg_copy(&st0, &tmp);
882 tmptag = st0_tag;
883 N = (expdif & 0x0000001f) + 32; /* This choice gives results
884 identical to an AMD 486 */
885 setexponent16(&tmp, N);
886 exp_1 = exponent16(&st1);
887 setexponent16(&st1, 0);
888 expdif -= N;
889
890 sign = getsign(&tmp) ^ st1_sign;
891 tag =
892 FPU_u_div(&tmp, &st1, &tmp,
893 PR_64_BITS | RC_CHOP | 0x3f, sign);
894 setsign(&tmp, sign);
895
896 FPU_round_to_int(&tmp, tag); /* Fortunately, this can't
897 overflow to 2^64 */
898
899 rem_kernel(significand(&st0),
900 &significand(&tmp),
901 significand(&st1),
902 significand(&tmp), exponent(&tmp)
903 );
904 setexponent16(&tmp, exp_1 + expdif);
905
906 /* It is possible for the operation to be complete here.
907 What does the IEEE standard say? The Intel 80486 manual
908 implies that the operation will never be completed at this
909 point, and the behaviour of a real 80486 confirms this.
910 */
911 if (!(tmp.sigh | tmp.sigl)) {
912 /* The result is zero */
913 control_word = old_cw;
914 partial_status = saved_status;
915 FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
916 setsign(&st0, st0_sign);
917#ifdef PECULIAR_486
918 setcc(SW_C2);
919#else
920 setcc(0);
921#endif /* PECULIAR_486 */
922 return;
923 }
924 cc = SW_C2;
925 }
926
927 control_word = old_cw;
928 partial_status = saved_status;
929 tag = FPU_normalize_nuo(&tmp);
930 reg_copy(&tmp, st0_ptr);
931
932 /* The only condition to be looked for is underflow,
933 and it can occur here only if underflow is unmasked. */
934 if ((exponent16(&tmp) <= EXP_UNDER) && (tag != TAG_Zero)
935 && !(control_word & CW_Underflow)) {
936 setcc(cc);
937 tag = arith_underflow(st0_ptr);
938 setsign(st0_ptr, st0_sign);
939 FPU_settag0(tag);
940 return;
941 } else if ((exponent16(&tmp) > EXP_UNDER) || (tag == TAG_Zero)) {
942 stdexp(st0_ptr);
943 setsign(st0_ptr, st0_sign);
944 } else {
945 tag =
946 FPU_round(st0_ptr, 0, 0, FULL_PRECISION, st0_sign);
947 }
948 FPU_settag0(tag);
949 setcc(cc);
950
951 return;
952 }
953
954 if (st0_tag == TAG_Special)
955 st0_tag = FPU_Special(st0_ptr);
956 if (st1_tag == TAG_Special)
957 st1_tag = FPU_Special(st1_ptr);
958
959 if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
960 || ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
961 || ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
962 if (denormal_operand() < 0)
963 return;
964 goto fprem_valid;
965 } else if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
966 FPU_stack_underflow();
967 return;
968 } else if (st0_tag == TAG_Zero) {
969 if (st1_tag == TAG_Valid) {
970 setcc(0);
971 return;
972 } else if (st1_tag == TW_Denormal) {
973 if (denormal_operand() < 0)
974 return;
975 setcc(0);
976 return;
977 } else if (st1_tag == TAG_Zero) {
978 arith_invalid(0);
979 return;
980 } /* fprem(?,0) always invalid */
981 else if (st1_tag == TW_Infinity) {
982 setcc(0);
983 return;
984 }
985 } else if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
986 if (st1_tag == TAG_Zero) {
987 arith_invalid(0); /* fprem(Valid,Zero) is invalid */
988 return;
989 } else if (st1_tag != TW_NaN) {
990 if (((st0_tag == TW_Denormal)
991 || (st1_tag == TW_Denormal))
992 && (denormal_operand() < 0))
993 return;
994
995 if (st1_tag == TW_Infinity) {
996 /* fprem(Valid,Infinity) is o.k. */
997 setcc(0);
998 return;
999 }
1000 }
1001 } else if (st0_tag == TW_Infinity) {
1002 if (st1_tag != TW_NaN) {
1003 arith_invalid(0); /* fprem(Infinity,?) is invalid */
1004 return;
1005 }
1006 }
1007
1008 /* One of the registers must contain a NaN if we got here. */
1009
1010#ifdef PARANOID
1011 if ((st0_tag != TW_NaN) && (st1_tag != TW_NaN))
1012 EXCEPTION(EX_INTERNAL | 0x118);
1013#endif /* PARANOID */
1014
1015 real_2op_NaN(st1_ptr, st1_tag, 0, st1_ptr);
1016
1017}
1018
1019/* ST(1) <- ST(1) * log ST; pop ST */
1020static void fyl2x(FPU_REG *st0_ptr, u_char st0_tag)
1021{
1022 FPU_REG *st1_ptr = &st(1), exponent;
1023 u_char st1_tag = FPU_gettagi(1);
1024 u_char sign;
1025 int e, tag;
1026
1027 clear_C1();
1028
1029 if ((st0_tag == TAG_Valid) && (st1_tag == TAG_Valid)) {
1030 both_valid:
1031 /* Both regs are Valid or Denormal */
1032 if (signpositive(st0_ptr)) {
1033 if (st0_tag == TW_Denormal)
1034 FPU_to_exp16(st0_ptr, st0_ptr);
1035 else
1036 /* Convert st(0) for internal use. */
1037 setexponent16(st0_ptr, exponent(st0_ptr));
1038
1039 if ((st0_ptr->sigh == 0x80000000)
1040 && (st0_ptr->sigl == 0)) {
1041 /* Special case. The result can be precise. */
1042 u_char esign;
1043 e = exponent16(st0_ptr);
1044 if (e >= 0) {
1045 exponent.sigh = e;
1046 esign = SIGN_POS;
1047 } else {
1048 exponent.sigh = -e;
1049 esign = SIGN_NEG;
1050 }
1051 exponent.sigl = 0;
1052 setexponent16(&exponent, 31);
1053 tag = FPU_normalize_nuo(&exponent);
1054 stdexp(&exponent);
1055 setsign(&exponent, esign);
1056 tag =
1057 FPU_mul(&exponent, tag, 1, FULL_PRECISION);
1058 if (tag >= 0)
1059 FPU_settagi(1, tag);
1060 } else {
1061 /* The usual case */
1062 sign = getsign(st1_ptr);
1063 if (st1_tag == TW_Denormal)
1064 FPU_to_exp16(st1_ptr, st1_ptr);
1065 else
1066 /* Convert st(1) for internal use. */
1067 setexponent16(st1_ptr,
1068 exponent(st1_ptr));
1069 poly_l2(st0_ptr, st1_ptr, sign);
1070 }
1071 } else {
1072 /* negative */
1073 if (arith_invalid(1) < 0)
1074 return;
1075 }
1076
1077 FPU_pop();
1078
1079 return;
1080 }
1081
1082 if (st0_tag == TAG_Special)
1083 st0_tag = FPU_Special(st0_ptr);
1084 if (st1_tag == TAG_Special)
1085 st1_tag = FPU_Special(st1_ptr);
1086
1087 if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
1088 FPU_stack_underflow_pop(1);
1089 return;
1090 } else if ((st0_tag <= TW_Denormal) && (st1_tag <= TW_Denormal)) {
1091 if (st0_tag == TAG_Zero) {
1092 if (st1_tag == TAG_Zero) {
1093 /* Both args zero is invalid */
1094 if (arith_invalid(1) < 0)
1095 return;
1096 } else {
1097 u_char sign;
1098 sign = getsign(st1_ptr) ^ SIGN_NEG;
1099 if (FPU_divide_by_zero(1, sign) < 0)
1100 return;
1101
1102 setsign(st1_ptr, sign);
1103 }
1104 } else if (st1_tag == TAG_Zero) {
1105 /* st(1) contains zero, st(0) valid <> 0 */
1106 /* Zero is the valid answer */
1107 sign = getsign(st1_ptr);
1108
1109 if (signnegative(st0_ptr)) {
1110 /* log(negative) */
1111 if (arith_invalid(1) < 0)
1112 return;
1113 } else if ((st0_tag == TW_Denormal)
1114 && (denormal_operand() < 0))
1115 return;
1116 else {
1117 if (exponent(st0_ptr) < 0)
1118 sign ^= SIGN_NEG;
1119
1120 FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
1121 setsign(st1_ptr, sign);
1122 }
1123 } else {
1124 /* One or both operands are denormals. */
1125 if (denormal_operand() < 0)
1126 return;
1127 goto both_valid;
1128 }
1129 } else if ((st0_tag == TW_NaN) || (st1_tag == TW_NaN)) {
1130 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1131 return;
1132 }
1133 /* One or both arg must be an infinity */
1134 else if (st0_tag == TW_Infinity) {
1135 if ((signnegative(st0_ptr)) || (st1_tag == TAG_Zero)) {
1136 /* log(-infinity) or 0*log(infinity) */
1137 if (arith_invalid(1) < 0)
1138 return;
1139 } else {
1140 u_char sign = getsign(st1_ptr);
1141
1142 if ((st1_tag == TW_Denormal)
1143 && (denormal_operand() < 0))
1144 return;
1145
1146 FPU_copy_to_reg1(&CONST_INF, TAG_Special);
1147 setsign(st1_ptr, sign);
1148 }
1149 }
1150 /* st(1) must be infinity here */
1151 else if (((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal))
1152 && (signpositive(st0_ptr))) {
1153 if (exponent(st0_ptr) >= 0) {
1154 if ((exponent(st0_ptr) == 0) &&
1155 (st0_ptr->sigh == 0x80000000) &&
1156 (st0_ptr->sigl == 0)) {
1157 /* st(0) holds 1.0 */
1158 /* infinity*log(1) */
1159 if (arith_invalid(1) < 0)
1160 return;
1161 }
1162 /* else st(0) is positive and > 1.0 */
1163 } else {
1164 /* st(0) is positive and < 1.0 */
1165
1166 if ((st0_tag == TW_Denormal)
1167 && (denormal_operand() < 0))
1168 return;
1169
1170 changesign(st1_ptr);
1171 }
1172 } else {
1173 /* st(0) must be zero or negative */
1174 if (st0_tag == TAG_Zero) {
1175 /* This should be invalid, but a real 80486 is happy with it. */
1176
1177#ifndef PECULIAR_486
1178 sign = getsign(st1_ptr);
1179 if (FPU_divide_by_zero(1, sign) < 0)
1180 return;
1181#endif /* PECULIAR_486 */
1182
1183 changesign(st1_ptr);
1184 } else if (arith_invalid(1) < 0) /* log(negative) */
1185 return;
1186 }
1187
1188 FPU_pop();
1189}
1190
1191static void fpatan(FPU_REG *st0_ptr, u_char st0_tag)
1192{
1193 FPU_REG *st1_ptr = &st(1);
1194 u_char st1_tag = FPU_gettagi(1);
1195 int tag;
1196
1197 clear_C1();
1198 if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1199 valid_atan:
1200
1201 poly_atan(st0_ptr, st0_tag, st1_ptr, st1_tag);
1202
1203 FPU_pop();
1204
1205 return;
1206 }
1207
1208 if (st0_tag == TAG_Special)
1209 st0_tag = FPU_Special(st0_ptr);
1210 if (st1_tag == TAG_Special)
1211 st1_tag = FPU_Special(st1_ptr);
1212
1213 if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
1214 || ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
1215 || ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
1216 if (denormal_operand() < 0)
1217 return;
1218
1219 goto valid_atan;
1220 } else if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
1221 FPU_stack_underflow_pop(1);
1222 return;
1223 } else if ((st0_tag == TW_NaN) || (st1_tag == TW_NaN)) {
1224 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) >= 0)
1225 FPU_pop();
1226 return;
1227 } else if ((st0_tag == TW_Infinity) || (st1_tag == TW_Infinity)) {
1228 u_char sign = getsign(st1_ptr);
1229 if (st0_tag == TW_Infinity) {
1230 if (st1_tag == TW_Infinity) {
1231 if (signpositive(st0_ptr)) {
1232 FPU_copy_to_reg1(&CONST_PI4, TAG_Valid);
1233 } else {
1234 setpositive(st1_ptr);
1235 tag =
1236 FPU_u_add(&CONST_PI4, &CONST_PI2,
1237 st1_ptr, FULL_PRECISION,
1238 SIGN_POS,
1239 exponent(&CONST_PI4),
1240 exponent(&CONST_PI2));
1241 if (tag >= 0)
1242 FPU_settagi(1, tag);
1243 }
1244 } else {
1245 if ((st1_tag == TW_Denormal)
1246 && (denormal_operand() < 0))
1247 return;
1248
1249 if (signpositive(st0_ptr)) {
1250 FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
1251 setsign(st1_ptr, sign); /* An 80486 preserves the sign */
1252 FPU_pop();
1253 return;
1254 } else {
1255 FPU_copy_to_reg1(&CONST_PI, TAG_Valid);
1256 }
1257 }
1258 } else {
1259 /* st(1) is infinity, st(0) not infinity */
1260 if ((st0_tag == TW_Denormal)
1261 && (denormal_operand() < 0))
1262 return;
1263
1264 FPU_copy_to_reg1(&CONST_PI2, TAG_Valid);
1265 }
1266 setsign(st1_ptr, sign);
1267 } else if (st1_tag == TAG_Zero) {
1268 /* st(0) must be valid or zero */
1269 u_char sign = getsign(st1_ptr);
1270
1271 if ((st0_tag == TW_Denormal) && (denormal_operand() < 0))
1272 return;
1273
1274 if (signpositive(st0_ptr)) {
1275 /* An 80486 preserves the sign */
1276 FPU_pop();
1277 return;
1278 }
1279
1280 FPU_copy_to_reg1(&CONST_PI, TAG_Valid);
1281 setsign(st1_ptr, sign);
1282 } else if (st0_tag == TAG_Zero) {
1283 /* st(1) must be TAG_Valid here */
1284 u_char sign = getsign(st1_ptr);
1285
1286 if ((st1_tag == TW_Denormal) && (denormal_operand() < 0))
1287 return;
1288
1289 FPU_copy_to_reg1(&CONST_PI2, TAG_Valid);
1290 setsign(st1_ptr, sign);
1291 }
1292#ifdef PARANOID
1293 else
1294 EXCEPTION(EX_INTERNAL | 0x125);
1295#endif /* PARANOID */
1296
1297 FPU_pop();
1298 set_precision_flag_up(); /* We do not really know if up or down */
1299}
1300
1301static void fprem(FPU_REG *st0_ptr, u_char st0_tag)
1302{
1303 do_fprem(st0_ptr, st0_tag, RC_CHOP);
1304}
1305
1306static void fprem1(FPU_REG *st0_ptr, u_char st0_tag)
1307{
1308 do_fprem(st0_ptr, st0_tag, RC_RND);
1309}
1310
1311static void fyl2xp1(FPU_REG *st0_ptr, u_char st0_tag)
1312{
1313 u_char sign, sign1;
1314 FPU_REG *st1_ptr = &st(1), a, b;
1315 u_char st1_tag = FPU_gettagi(1);
1316
1317 clear_C1();
1318 if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1319 valid_yl2xp1:
1320
1321 sign = getsign(st0_ptr);
1322 sign1 = getsign(st1_ptr);
1323
1324 FPU_to_exp16(st0_ptr, &a);
1325 FPU_to_exp16(st1_ptr, &b);
1326
1327 if (poly_l2p1(sign, sign1, &a, &b, st1_ptr))
1328 return;
1329
1330 FPU_pop();
1331 return;
1332 }
1333
1334 if (st0_tag == TAG_Special)
1335 st0_tag = FPU_Special(st0_ptr);
1336 if (st1_tag == TAG_Special)
1337 st1_tag = FPU_Special(st1_ptr);
1338
1339 if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
1340 || ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
1341 || ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
1342 if (denormal_operand() < 0)
1343 return;
1344
1345 goto valid_yl2xp1;
1346 } else if ((st0_tag == TAG_Empty) | (st1_tag == TAG_Empty)) {
1347 FPU_stack_underflow_pop(1);
1348 return;
1349 } else if (st0_tag == TAG_Zero) {
1350 switch (st1_tag) {
1351 case TW_Denormal:
1352 if (denormal_operand() < 0)
1353 return;
1354
1355 case TAG_Zero:
1356 case TAG_Valid:
1357 setsign(st0_ptr, getsign(st0_ptr) ^ getsign(st1_ptr));
1358 FPU_copy_to_reg1(st0_ptr, st0_tag);
1359 break;
1360
1361 case TW_Infinity:
1362 /* Infinity*log(1) */
1363 if (arith_invalid(1) < 0)
1364 return;
1365 break;
1366
1367 case TW_NaN:
1368 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1369 return;
1370 break;
1371
1372 default:
1373#ifdef PARANOID
1374 EXCEPTION(EX_INTERNAL | 0x116);
1375 return;
1376#endif /* PARANOID */
1377 break;
1378 }
1379 } else if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
1380 switch (st1_tag) {
1381 case TAG_Zero:
1382 if (signnegative(st0_ptr)) {
1383 if (exponent(st0_ptr) >= 0) {
1384 /* st(0) holds <= -1.0 */
1385#ifdef PECULIAR_486 /* Stupid 80486 doesn't worry about log(negative). */
1386 changesign(st1_ptr);
1387#else
1388 if (arith_invalid(1) < 0)
1389 return;
1390#endif /* PECULIAR_486 */
1391 } else if ((st0_tag == TW_Denormal)
1392 && (denormal_operand() < 0))
1393 return;
1394 else
1395 changesign(st1_ptr);
1396 } else if ((st0_tag == TW_Denormal)
1397 && (denormal_operand() < 0))
1398 return;
1399 break;
1400
1401 case TW_Infinity:
1402 if (signnegative(st0_ptr)) {
1403 if ((exponent(st0_ptr) >= 0) &&
1404 !((st0_ptr->sigh == 0x80000000) &&
1405 (st0_ptr->sigl == 0))) {
1406 /* st(0) holds < -1.0 */
1407#ifdef PECULIAR_486 /* Stupid 80486 doesn't worry about log(negative). */
1408 changesign(st1_ptr);
1409#else
1410 if (arith_invalid(1) < 0)
1411 return;
1412#endif /* PECULIAR_486 */
1413 } else if ((st0_tag == TW_Denormal)
1414 && (denormal_operand() < 0))
1415 return;
1416 else
1417 changesign(st1_ptr);
1418 } else if ((st0_tag == TW_Denormal)
1419 && (denormal_operand() < 0))
1420 return;
1421 break;
1422
1423 case TW_NaN:
1424 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1425 return;
1426 }
1427
1428 } else if (st0_tag == TW_NaN) {
1429 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1430 return;
1431 } else if (st0_tag == TW_Infinity) {
1432 if (st1_tag == TW_NaN) {
1433 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1434 return;
1435 } else if (signnegative(st0_ptr)) {
1436#ifndef PECULIAR_486
1437 /* This should have higher priority than denormals, but... */
1438 if (arith_invalid(1) < 0) /* log(-infinity) */
1439 return;
1440#endif /* PECULIAR_486 */
1441 if ((st1_tag == TW_Denormal)
1442 && (denormal_operand() < 0))
1443 return;
1444#ifdef PECULIAR_486
1445 /* Denormal operands actually get higher priority */
1446 if (arith_invalid(1) < 0) /* log(-infinity) */
1447 return;
1448#endif /* PECULIAR_486 */
1449 } else if (st1_tag == TAG_Zero) {
1450 /* log(infinity) */
1451 if (arith_invalid(1) < 0)
1452 return;
1453 }
1454
1455 /* st(1) must be valid here. */
1456
1457 else if ((st1_tag == TW_Denormal) && (denormal_operand() < 0))
1458 return;
1459
1460 /* The Manual says that log(Infinity) is invalid, but a real
1461 80486 sensibly says that it is o.k. */
1462 else {
1463 u_char sign = getsign(st1_ptr);
1464 FPU_copy_to_reg1(&CONST_INF, TAG_Special);
1465 setsign(st1_ptr, sign);
1466 }
1467 }
1468#ifdef PARANOID
1469 else {
1470 EXCEPTION(EX_INTERNAL | 0x117);
1471 return;
1472 }
1473#endif /* PARANOID */
1474
1475 FPU_pop();
1476 return;
1477
1478}
1479
1480static void fscale(FPU_REG *st0_ptr, u_char st0_tag)
1481{
1482 FPU_REG *st1_ptr = &st(1);
1483 u_char st1_tag = FPU_gettagi(1);
1484 int old_cw = control_word;
1485 u_char sign = getsign(st0_ptr);
1486
1487 clear_C1();
1488 if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1489 long scale;
1490 FPU_REG tmp;
1491
1492 /* Convert register for internal use. */
1493 setexponent16(st0_ptr, exponent(st0_ptr));
1494
1495 valid_scale:
1496
1497 if (exponent(st1_ptr) > 30) {
1498 /* 2^31 is far too large, would require 2^(2^30) or 2^(-2^30) */
1499
1500 if (signpositive(st1_ptr)) {
1501 EXCEPTION(EX_Overflow);
1502 FPU_copy_to_reg0(&CONST_INF, TAG_Special);
1503 } else {
1504 EXCEPTION(EX_Underflow);
1505 FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
1506 }
1507 setsign(st0_ptr, sign);
1508 return;
1509 }
1510
1511 control_word &= ~CW_RC;
1512 control_word |= RC_CHOP;
1513 reg_copy(st1_ptr, &tmp);
1514 FPU_round_to_int(&tmp, st1_tag); /* This can never overflow here */
1515 control_word = old_cw;
1516 scale = signnegative(st1_ptr) ? -tmp.sigl : tmp.sigl;
1517 scale += exponent16(st0_ptr);
1518
1519 setexponent16(st0_ptr, scale);
1520
1521 /* Use FPU_round() to properly detect under/overflow etc */
1522 FPU_round(st0_ptr, 0, 0, control_word, sign);
1523
1524 return;
1525 }
1526
1527 if (st0_tag == TAG_Special)
1528 st0_tag = FPU_Special(st0_ptr);
1529 if (st1_tag == TAG_Special)
1530 st1_tag = FPU_Special(st1_ptr);
1531
1532 if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
1533 switch (st1_tag) {
1534 case TAG_Valid:
1535 /* st(0) must be a denormal */
1536 if ((st0_tag == TW_Denormal)
1537 && (denormal_operand() < 0))
1538 return;
1539
1540 FPU_to_exp16(st0_ptr, st0_ptr); /* Will not be left on stack */
1541 goto valid_scale;
1542
1543 case TAG_Zero:
1544 if (st0_tag == TW_Denormal)
1545 denormal_operand();
1546 return;
1547
1548 case TW_Denormal:
1549 denormal_operand();
1550 return;
1551
1552 case TW_Infinity:
1553 if ((st0_tag == TW_Denormal)
1554 && (denormal_operand() < 0))
1555 return;
1556
1557 if (signpositive(st1_ptr))
1558 FPU_copy_to_reg0(&CONST_INF, TAG_Special);
1559 else
1560 FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
1561 setsign(st0_ptr, sign);
1562 return;
1563
1564 case TW_NaN:
1565 real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1566 return;
1567 }
1568 } else if (st0_tag == TAG_Zero) {
1569 switch (st1_tag) {
1570 case TAG_Valid:
1571 case TAG_Zero:
1572 return;
1573
1574 case TW_Denormal:
1575 denormal_operand();
1576 return;
1577
1578 case TW_Infinity:
1579 if (signpositive(st1_ptr))
1580 arith_invalid(0); /* Zero scaled by +Infinity */
1581 return;
1582
1583 case TW_NaN:
1584 real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1585 return;
1586 }
1587 } else if (st0_tag == TW_Infinity) {
1588 switch (st1_tag) {
1589 case TAG_Valid:
1590 case TAG_Zero:
1591 return;
1592
1593 case TW_Denormal:
1594 denormal_operand();
1595 return;
1596
1597 case TW_Infinity:
1598 if (signnegative(st1_ptr))
1599 arith_invalid(0); /* Infinity scaled by -Infinity */
1600 return;
1601
1602 case TW_NaN:
1603 real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1604 return;
1605 }
1606 } else if (st0_tag == TW_NaN) {
1607 if (st1_tag != TAG_Empty) {
1608 real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1609 return;
1610 }
1611 }
1612#ifdef PARANOID
1613 if (!((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty))) {
1614 EXCEPTION(EX_INTERNAL | 0x115);
1615 return;
1616 }
1617#endif
1618
1619 /* At least one of st(0), st(1) must be empty */
1620 FPU_stack_underflow();
1621
1622}
1623
1624/*---------------------------------------------------------------------------*/
1625
1626static FUNC_ST0 const trig_table_a[] = {
1627 f2xm1, fyl2x, fptan, fpatan,
1628 fxtract, fprem1, (FUNC_ST0) fdecstp, (FUNC_ST0) fincstp
1629};
1630
1631void FPU_triga(void)
1632{
1633 (trig_table_a[FPU_rm]) (&st(0), FPU_gettag0());
1634}
1635
1636static FUNC_ST0 const trig_table_b[] = {
1637 fprem, fyl2xp1, fsqrt_, fsincos, frndint_, fscale, (FUNC_ST0) fsin, fcos
1638};
1639
1640void FPU_trigb(void)
1641{
1642 (trig_table_b[FPU_rm]) (&st(0), FPU_gettag0());
1643}