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

DEFINITIONS
