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

DEFINITIONS
