kernel/FPU-emu/fpu

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

DEFINITIONS
