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

DEFINITIONS
