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