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

DEFINITIONS
