arch/i386/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,
  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 }