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