arch/x86/math-emu/fpu_trig.c

   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, 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 */
  33 static 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 */
 151 static 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
 179 static 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
 191 static 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
 232 static 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
 275 static 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
 354 static 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
 435 static void fdecstp(void)
 436 {
 437         clear_C1();
 438         top--;
 439 }
 440
 441 static void fincstp(void)
 442 {
 443         clear_C1();
 444         top++;
 445 }
 446
 447 static 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
 505 static 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
 549 static 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
 610 static 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
 680 static void fcos(FPU_REG *st0_ptr, u_char st0_tag)
 681 {
 682         f_cos(st0_ptr, st0_tag);
 683 }
 684
 685 static 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  */
 745 static 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. */
 778 static 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 */
1020 static 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
1191 static 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
1301 static void fprem(FPU_REG *st0_ptr, u_char st0_tag)
1302 {
1303         do_fprem(st0_ptr, st0_tag, RC_CHOP);
1304 }
1305
1306 static void fprem1(FPU_REG *st0_ptr, u_char st0_tag)
1307 {
1308         do_fprem(st0_ptr, st0_tag, RC_RND);
1309 }
1310
1311 static 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
1480 static 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
1626 static FUNC_ST0 const trig_table_a[] = {
1627         f2xm1, fyl2x, fptan, fpatan,
1628         fxtract, fprem1, (FUNC_ST0) fdecstp, (FUNC_ST0) fincstp
1629 };
1630
1631 void FPU_triga(void)
1632 {
1633         (trig_table_a[FPU_rm]) (&st(0), FPU_gettag0());
1634 }
1635
1636 static FUNC_ST0 const trig_table_b[] = {
1637         fprem, fyl2xp1, fsqrt_, fsincos, frndint_, fscale, (FUNC_ST0) fsin, fcos
1638 };
1639
1640 void FPU_trigb(void)
1641 {
1642         (trig_table_b[FPU_rm]) (&st(0), FPU_gettag0());
1643 }