kernel/FPU-emu/fpu

/* */
This source file includes following definitions.
trig_arg
convert_l2reg
single_arg_error
f2xm1
fptan
fxtract
fdecstp
fincstp
fsqrt_
frndint_
fsin
f_cos
fcos
fsincos
fprem_kernel
fyl2x
fpatan
fprem
fprem1
fyl2xp1
fscale
trig_a
trig_b
   1 /*---------------------------------------------------------------------------+
   2  |  fpu_trig.c                                                               |
   3  |                                                                           |
   4  | Implementation of the FPU "transcendental" functions.                     |
   5  |                                                                           |
   6  | Copyright (C) 1992,1993                                                   |
   7  |                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |
   8  |                       Australia.  E-mail apm233m@vaxc.cc.monash.edu.au    |
   9  |                                                                           |
  10  |                                                                           |
  11  +---------------------------------------------------------------------------*/
  12 
  13 #include "fpu_system.h"
  14 #include "exception.h"
  15 #include "fpu_emu.h"
  16 #include "status_w.h"
  17 #include "control_w.h"
  18 #include "reg_constant.h"       
  19 
  20 
  21 
  22 static int trig_arg(FPU_REG *X)
     /*  */
  23 {
  24   FPU_REG tmp, quot;
  25   int rv;
  26   long long q;
  27   int old_cw = control_word;
  28 
  29   control_word &= ~CW_RC;
  30   control_word |= RC_CHOP;
  31   
  32   reg_move(X, &quot);
  33   reg_div(&quot, &CONST_PI2, &quot, FULL_PRECISION);
  34 
  35   reg_move(&quot, &tmp);
  36   round_to_int(&tmp);
  37   if ( tmp.sigh & 0x80000000 )
  38     return -1;              /* |Arg| is >= 2^63 */
  39   tmp.exp = EXP_BIAS + 63;
  40   q = *(long long *)&(tmp.sigl);
  41   normalize(&tmp);
  42 
  43   reg_sub(&quot, &tmp, X, FULL_PRECISION);
  44   rv = q & 7;
  45   
  46   control_word = old_cw;
  47   return rv;;
  48 }
  49 
  50 
  51 /* Convert a long to register */
  52 void convert_l2reg(long *arg, FPU_REG *dest)
     /*  */
  53 {
  54   long num = *arg;
  55   
  56   if (num == 0)
  57     { reg_move(&CONST_Z, dest); return; }
  58 
  59   if (num > 0)
  60     dest->sign = SIGN_POS;
  61   else
  62     { num = -num; dest->sign = SIGN_NEG; }
  63   
  64   dest->sigh = num;
  65   dest->sigl = 0;
  66   dest->exp = EXP_BIAS + 31;
  67   dest->tag = TW_Valid;
  68   normalize(dest);
  69 }
  70 
  71 
  72 static void single_arg_error(void)
     /*  */
  73 {
  74   switch ( FPU_st0_tag )
  75     {
  76     case TW_NaN:
  77       if ( !(FPU_st0_ptr->sigh & 0x40000000) )   /* Signaling ? */
  78         {
  79           EXCEPTION(EX_Invalid);
  80           /* Convert to a QNaN */
  81           FPU_st0_ptr->sigh |= 0x40000000;
  82         }
  83       break;              /* return with a NaN in st(0) */
  84     case TW_Empty:
  85       stack_underflow();  /* Puts a QNaN in st(0) */
  86       break;
  87 #ifdef PARANOID
  88     default:
  89       EXCEPTION(EX_INTERNAL|0x0112);
  90 #endif PARANOID
  91     }
  92 }
  93 
  94 
  95 /*---------------------------------------------------------------------------*/
  96 
  97 static void f2xm1(void)
     /*  */
  98 {
  99   switch ( FPU_st0_tag )
 100     {
 101     case TW_Valid:
 102       {
 103         FPU_REG rv, tmp;
 104 
 105 #ifdef DENORM_OPERAND
 106       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
 107         return;
 108 #endif DENORM_OPERAND
 109         
 110         if ( FPU_st0_ptr->sign == SIGN_POS )
 111           {
 112             /* poly_2xm1(x) requires 0 < x < 1. */
 113             if ( poly_2xm1(FPU_st0_ptr, &rv) )
 114               return;                  /* error */
 115             reg_mul(&rv, FPU_st0_ptr, FPU_st0_ptr, FULL_PRECISION);
 116           }
 117         else
 118           {
 119 /* **** Should change poly_2xm1() to at least handle numbers near 0 */
 120             /* poly_2xm1(x) doesn't handle negative numbers. */
 121             /* So we compute (poly_2xm1(x+1)-1)/2, for -1 < x < 0 */
 122             reg_add(FPU_st0_ptr, &CONST_1, &tmp, FULL_PRECISION);
 123             poly_2xm1(&tmp, &rv);
 124             reg_mul(&rv, &tmp, &tmp, FULL_PRECISION);
 125             reg_sub(&tmp, &CONST_1, FPU_st0_ptr, FULL_PRECISION);
 126             FPU_st0_ptr->exp--;
 127             if ( FPU_st0_ptr->exp <= EXP_UNDER )
 128               arith_underflow(FPU_st0_ptr);
 129           }
 130         return;
 131       }
 132     case TW_Zero:
 133       return;
 134     case TW_Infinity:
 135       if ( FPU_st0_ptr->sign == SIGN_NEG )
 136         {
 137           /* -infinity gives -1 (p16-10) */
 138           reg_move(&CONST_1, FPU_st0_ptr);
 139           FPU_st0_ptr->sign = SIGN_NEG;
 140         }
 141       return;
 142     default:
 143       single_arg_error();
 144     }
 145 }
 146 
 147 static void fptan(void)
     /*  */
 148 {
 149   FPU_REG *st_new_ptr;
 150   int q;
 151   char arg_sign = FPU_st0_ptr->sign;
 152 
 153   if ( STACK_OVERFLOW )
 154     { stack_overflow(); return; }
 155 
 156   switch ( FPU_st0_tag )
 157     {
 158     case TW_Valid:
 159 
 160 #ifdef DENORM_OPERAND
 161       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
 162         return;
 163 #endif DENORM_OPERAND
 164 
 165       FPU_st0_ptr->sign = SIGN_POS;
 166       if ( (q = trig_arg(FPU_st0_ptr)) != -1 )
 167         {
 168           if (q & 1)
 169             reg_sub(&CONST_1, FPU_st0_ptr, FPU_st0_ptr, FULL_PRECISION);
 170           
 171           poly_tan(FPU_st0_ptr, FPU_st0_ptr);
 172 
 173           FPU_st0_ptr->sign = (q & 1) ^ arg_sign;
 174           
 175           if ( FPU_st0_ptr->exp <= EXP_UNDER )
 176             arith_underflow(FPU_st0_ptr);
 177           
 178           push();
 179           reg_move(&CONST_1, FPU_st0_ptr);
 180           setcc(0);
 181         }
 182       else
 183         {
 184           /* Operand is out of range */
 185           setcc(SW_C2);
 186           FPU_st0_ptr->sign = arg_sign;         /* restore st(0) */
 187           return;
 188         }
 189       break;
 190     case TW_Infinity:
 191       /* Operand is out of range */
 192       setcc(SW_C2);
 193       FPU_st0_ptr->sign = arg_sign;         /* restore st(0) */
 194       return;
 195     case TW_Zero:
 196       push();
 197       reg_move(&CONST_1, FPU_st0_ptr);
 198       setcc(0);
 199       break;
 200     default:
 201       single_arg_error();
 202       break;
 203     }
 204 }
 205 
 206 
 207 static void fxtract(void)
     /*  */
 208 {
 209   FPU_REG *st_new_ptr;
 210   register FPU_REG *st1_ptr = FPU_st0_ptr;  /* anticipate */
 211 
 212   if ( STACK_OVERFLOW )
 213     {  stack_overflow(); return; }
 214 
 215   if ( !(FPU_st0_tag ^ TW_Valid) )
 216     {
 217       long e;
 218 
 219 #ifdef DENORM_OPERAND
 220       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
 221         return;
 222 #endif DENORM_OPERAND
 223           
 224       push();
 225       reg_move(st1_ptr, FPU_st0_ptr);
 226       FPU_st0_ptr->exp = EXP_BIAS;
 227       e = st1_ptr->exp - EXP_BIAS;
 228       convert_l2reg(&e, st1_ptr);
 229       return;
 230     }
 231   else if ( FPU_st0_tag == TW_Zero )
 232     {
 233       char sign = FPU_st0_ptr->sign;
 234       divide_by_zero(SIGN_NEG, FPU_st0_ptr);
 235       push();
 236       reg_move(&CONST_Z, FPU_st0_ptr);
 237       FPU_st0_ptr->sign = sign;
 238       return;
 239     }
 240   else if ( FPU_st0_tag == TW_Infinity )
 241     {
 242       char sign = FPU_st0_ptr->sign;
 243       FPU_st0_ptr->sign = SIGN_POS;
 244       push();
 245       reg_move(&CONST_INF, FPU_st0_ptr);
 246       FPU_st0_ptr->sign = sign;
 247       return;
 248     }
 249   else if ( FPU_st0_tag == TW_NaN )
 250     {
 251       if ( !(FPU_st0_ptr->sigh & 0x40000000) )   /* Signaling ? */
 252         {
 253           EXCEPTION(EX_Invalid);
 254           /* Convert to a QNaN */
 255           FPU_st0_ptr->sigh |= 0x40000000;
 256         }
 257       push();
 258       reg_move(st1_ptr, FPU_st0_ptr);
 259       return;
 260     }
 261   else if ( FPU_st0_tag == TW_Empty )
 262     {
 263       /* Is this the correct behaviour? */
 264       if ( control_word & EX_Invalid )
 265         {
 266           stack_underflow();
 267           push();
 268           stack_underflow();
 269         }
 270       else
 271         EXCEPTION(EX_StackUnder);
 272     }
 273 #ifdef PARANOID
 274   else
 275     EXCEPTION(EX_INTERNAL | 0x119);
 276 #endif PARANOID
 277 }
 278 
 279 
 280 static void fdecstp(void)
     /*  */
 281 {
 282   top--;  /* FPU_st0_ptr will be fixed in math_emulate() before the next instr */
 283 }
 284 
 285 static void fincstp(void)
     /*  */
 286 {
 287   top++;  /* FPU_st0_ptr will be fixed in math_emulate() before the next instr */
 288 }
 289 
 290 
 291 static void fsqrt_(void)
     /*  */
 292 {
 293   if ( !(FPU_st0_tag ^ TW_Valid) )
 294     {
 295       int expon;
 296       
 297       if (FPU_st0_ptr->sign == SIGN_NEG)
 298         {
 299           arith_invalid(FPU_st0_ptr);  /* sqrt(negative) is invalid */
 300           return;
 301         }
 302 
 303 #ifdef DENORM_OPERAND
 304       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
 305         return;
 306 #endif DENORM_OPERAND
 307 
 308       expon = FPU_st0_ptr->exp - EXP_BIAS;
 309       FPU_st0_ptr->exp = EXP_BIAS + (expon & 1);  /* make st(0) in  [1.0 .. 4.0) */
 310       
 311       wm_sqrt(FPU_st0_ptr, control_word);       /* Do the computation */
 312       
 313       FPU_st0_ptr->exp += expon >> 1;
 314       FPU_st0_ptr->sign = SIGN_POS;
 315     }
 316   else if ( FPU_st0_tag == TW_Zero )
 317     return;
 318   else if ( FPU_st0_tag == TW_Infinity )
 319     {
 320       if ( FPU_st0_ptr->sign == SIGN_NEG )
 321         arith_invalid(FPU_st0_ptr);  /* sqrt(-Infinity) is invalid */
 322       return;
 323     }
 324   else
 325     { single_arg_error(); return; }
 326 
 327 }
 328 
 329 
 330 static void frndint_(void)
     /*  */
 331 {
 332   if ( !(FPU_st0_tag ^ TW_Valid) )
 333     {
 334       if (FPU_st0_ptr->exp > EXP_BIAS+63)
 335         return;
 336 
 337 #ifdef DENORM_OPERAND
 338       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
 339         return;
 340 #endif DENORM_OPERAND
 341 
 342       round_to_int(FPU_st0_ptr);  /* Fortunately, this can't overflow to 2^64 */
 343       FPU_st0_ptr->exp = EXP_BIAS + 63;
 344       normalize(FPU_st0_ptr);
 345       return;
 346     }
 347   else if ( (FPU_st0_tag == TW_Zero) || (FPU_st0_tag == TW_Infinity) )
 348     return;
 349   else
 350     single_arg_error();
 351 }
 352 
 353 
 354 static void fsin(void)
     /*  */
 355 {
 356   char arg_sign = FPU_st0_ptr->sign;
 357 
 358   if ( FPU_st0_tag == TW_Valid )
 359     {
 360       int q;
 361       FPU_st0_ptr->sign = SIGN_POS;
 362 
 363 #ifdef DENORM_OPERAND
 364       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
 365         return;
 366 #endif DENORM_OPERAND
 367 
 368       if ( (q = trig_arg(FPU_st0_ptr)) != -1 )
 369         {
 370           FPU_REG rv;
 371           
 372           if (q & 1)
 373             reg_sub(&CONST_1, FPU_st0_ptr, FPU_st0_ptr, FULL_PRECISION);
 374           
 375           poly_sine(FPU_st0_ptr, &rv);
 376           
 377           setcc(0);
 378           if (q & 2)
 379             rv.sign ^= SIGN_POS ^ SIGN_NEG;
 380           rv.sign ^= arg_sign;
 381           reg_move(&rv, FPU_st0_ptr);
 382 
 383           if ( FPU_st0_ptr->exp <= EXP_UNDER )
 384             arith_underflow(FPU_st0_ptr);
 385 
 386           set_precision_flag_up();  /* We do not really know if up or down */
 387 
 388           return;
 389         }
 390       else
 391         {
 392           /* Operand is out of range */
 393           setcc(SW_C2);
 394           FPU_st0_ptr->sign = arg_sign;         /* restore st(0) */
 395           return;
 396         }
 397     }
 398   else if ( FPU_st0_tag == TW_Zero )
 399     {
 400       setcc(0);
 401       return;
 402     }
 403   else if ( FPU_st0_tag == TW_Infinity )
 404     {
 405       /* Operand is out of range */
 406       setcc(SW_C2);
 407       FPU_st0_ptr->sign = arg_sign;         /* restore st(0) */
 408       return;
 409     }
 410   else
 411     single_arg_error();
 412 }
 413 
 414 
 415 static int f_cos(FPU_REG *arg)
     /*  */
 416 {
 417   char arg_sign = arg->sign;
 418 
 419   if ( arg->tag == TW_Valid )
 420     {
 421       int q;
 422 
 423 #ifdef DENORM_OPERAND
 424       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
 425         return 1;
 426 #endif DENORM_OPERAND
 427 
 428       arg->sign = SIGN_POS;
 429       if ( (q = trig_arg(arg)) != -1 )
 430         {
 431           FPU_REG rv;
 432           
 433           if ( !(q & 1) )
 434             reg_sub(&CONST_1, arg, arg, FULL_PRECISION);
 435           
 436           poly_sine(arg, &rv);
 437           
 438           setcc(0);
 439           if ((q+1) & 2)
 440             rv.sign ^= SIGN_POS ^ SIGN_NEG;
 441           reg_move(&rv, arg);
 442 
 443           set_precision_flag_up();  /* We do not really know if up or down */
 444           
 445           return 0;
 446         }
 447       else
 448         {
 449           /* Operand is out of range */
 450           setcc(SW_C2);
 451           arg->sign = arg_sign;         /* restore st(0) */
 452           return 1;
 453         }
 454     }
 455   else if ( arg->tag == TW_Zero )
 456     {
 457       reg_move(&CONST_1, arg);
 458       setcc(0);
 459       return 0;
 460     }
 461   else if ( FPU_st0_tag == TW_Infinity )
 462     {
 463       /* Operand is out of range */
 464       setcc(SW_C2);
 465       arg->sign = arg_sign;         /* restore st(0) */
 466       return 1;
 467     }
 468   else
 469     {
 470       single_arg_error();  /* requires arg == &st(0) */
 471       return 1;
 472     }
 473 }
 474 
 475 
 476 static void fcos(void)
     /*  */
 477 {
 478   f_cos(FPU_st0_ptr);
 479 }
 480 
 481 
 482 static void fsincos(void)
     /*  */
 483 {
 484   FPU_REG *st_new_ptr;
 485   FPU_REG arg;
 486           
 487   if ( STACK_OVERFLOW )
 488     { stack_overflow(); return; }
 489 
 490   reg_move(FPU_st0_ptr,&arg);
 491   if ( !f_cos(&arg) )
 492     {
 493       fsin();
 494       push();
 495       reg_move(&arg,FPU_st0_ptr);
 496     }
 497 
 498 }
 499 
 500 
 501 /*---------------------------------------------------------------------------*/
 502 /* The following all require two arguments: st(0) and st(1) */
 503 
 504 /* remainder of st(0) / st(1) */
 505 /* Assumes that st(0) and st(1) are both TW_Valid */
 506 static void fprem_kernel(int round)
     /*  */
 507 {
 508   FPU_REG *st1_ptr = &st(1);
 509   char st1_tag = st1_ptr->tag;
 510 
 511   if ( !((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid)) )
 512     {
 513       FPU_REG tmp;
 514       int old_cw = control_word;
 515       int expdif = FPU_st0_ptr->exp - (st1_ptr)->exp;
 516 
 517 #ifdef DENORM_OPERAND
 518           if ( ((FPU_st0_ptr->exp <= EXP_UNDER) ||
 519                 (st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()) )
 520             return;
 521 #endif DENORM_OPERAND
 522       
 523       control_word &= ~CW_RC;
 524       control_word |= round;
 525       
 526       if (expdif < 64)
 527         {
 528           /* This should be the most common case */
 529           long long q;
 530           int c = 0;
 531 
 532           reg_div(FPU_st0_ptr, st1_ptr, &tmp, FULL_PRECISION);
 533           
 534           round_to_int(&tmp);  /* Fortunately, this can't overflow to 2^64 */
 535           tmp.exp = EXP_BIAS + 63;
 536           q = *(long long *)&(tmp.sigl);
 537           normalize(&tmp);
 538           
 539           reg_mul(st1_ptr, &tmp, &tmp, FULL_PRECISION);
 540           reg_sub(FPU_st0_ptr, &tmp, FPU_st0_ptr, FULL_PRECISION);
 541           
 542           if (q&4) c |= SW_C3;
 543           if (q&2) c |= SW_C1;
 544           if (q&1) c |= SW_C0;
 545           
 546           setcc(c);
 547         }
 548       else
 549         {
 550           /* There is a large exponent difference ( >= 64 ) */
 551           int N_exp;
 552           
 553           reg_div(FPU_st0_ptr, st1_ptr, &tmp, FULL_PRECISION);
 554           /* N is 'a number between 32 and 63' (p26-113) */
 555           N_exp = (tmp.exp & 31) + 32;
 556           tmp.exp = EXP_BIAS + N_exp;
 557           
 558           round_to_int(&tmp);  /* Fortunately, this can't overflow to 2^64 */
 559           tmp.exp = EXP_BIAS + 63;
 560           normalize(&tmp);
 561           
 562           tmp.exp = EXP_BIAS + expdif - N_exp;
 563           
 564           reg_mul(st1_ptr, &tmp, &tmp, FULL_PRECISION);
 565           reg_sub(FPU_st0_ptr, &tmp, FPU_st0_ptr, FULL_PRECISION);
 566           
 567           setcc(SW_C2);
 568         }
 569       control_word = old_cw;
 570 
 571       if ( FPU_st0_ptr->exp <= EXP_UNDER )
 572         arith_underflow(FPU_st0_ptr);
 573       return;
 574     }
 575   else if ( (FPU_st0_tag == TW_Empty) | (st1_tag == TW_Empty) )
 576     { stack_underflow(); return; }
 577   else if ( FPU_st0_tag == TW_Zero )
 578     {
 579       if ( st1_tag == TW_Valid )
 580         {
 581 
 582 #ifdef DENORM_OPERAND
 583           if ( (st1_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
 584             return;
 585 #endif DENORM_OPERAND
 586 
 587           setcc(0); return;
 588         }
 589       else if ( st1_tag == TW_Zero )
 590         { arith_invalid(FPU_st0_ptr); return; } /* fprem(?,0) always invalid */
 591       else if ( st1_tag == TW_Infinity )
 592         { setcc(0); return; }
 593     }
 594   else if ( FPU_st0_tag == TW_Valid )
 595     {
 596       if ( st1_tag == TW_Zero )
 597         {
 598           arith_invalid(FPU_st0_ptr); /* fprem(Valid,Zero) is invalid */
 599           return;
 600         }
 601       else if ( st1_tag != TW_NaN )
 602         {
 603 #ifdef DENORM_OPERAND
 604           if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
 605             return;
 606 #endif DENORM_OPERAND
 607 
 608           if ( st1_tag == TW_Infinity )
 609             {
 610               /* fprem(Valid,Infinity) is o.k. */
 611               setcc(0); return;
 612             }
 613         }
 614     }
 615   else if ( FPU_st0_tag == TW_Infinity )
 616     {
 617       if ( st1_tag != TW_NaN )
 618         {
 619           arith_invalid(FPU_st0_ptr); /* fprem(Infinity,?) is invalid */
 620           return;
 621         }
 622     }
 623 
 624   /* One of the registers must contain a NaN is we got here. */
 625 
 626 #ifdef PARANOID
 627   if ( (FPU_st0_tag != TW_NaN) && (st1_tag != TW_NaN) )
 628       EXCEPTION(EX_INTERNAL | 0x118);
 629 #endif PARANOID
 630 
 631   real_2op_NaN(FPU_st0_ptr, st1_ptr, FPU_st0_ptr);
 632 
 633 }
 634 
 635 
 636 /* ST(1) <- ST(1) * log ST;  pop ST */
 637 static void fyl2x(void)
     /*  */
 638 {
 639   FPU_REG *st1_ptr = &st(1);
 640   char st1_tag = st1_ptr->tag;
 641 
 642   if ( !((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid)) )
 643     {
 644       if ( FPU_st0_ptr->sign == SIGN_POS )
 645         {
 646           int saved_control, saved_status;
 647 
 648 #ifdef DENORM_OPERAND
 649           if ( ((FPU_st0_ptr->exp <= EXP_UNDER) ||
 650                 (st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()) )
 651             return;
 652 #endif DENORM_OPERAND
 653 
 654           /* We use the general purpose arithmetic,
 655              so we need to save these. */
 656           saved_status = status_word;
 657           saved_control = control_word;
 658           control_word = FULL_PRECISION;
 659 
 660           poly_l2(FPU_st0_ptr, FPU_st0_ptr);
 661 
 662           /* Enough of the basic arithmetic is done now */
 663           control_word = saved_control;
 664           status_word = saved_status;
 665 
 666           /* Let the multiply set the flags */
 667           reg_mul(FPU_st0_ptr, st1_ptr, st1_ptr, FULL_PRECISION);
 668 
 669           pop(); FPU_st0_ptr = &st(0);
 670         }
 671       else
 672         {
 673           /* negative   */
 674           pop(); FPU_st0_ptr = &st(0);
 675           arith_invalid(FPU_st0_ptr);  /* st(0) cannot be negative */
 676           return;
 677         }
 678     }
 679   else if ( (FPU_st0_tag == TW_Empty) || (st1_tag == TW_Empty) )
 680     {
 681       stack_underflow_pop(1);
 682       return;
 683     }
 684   else if ( (FPU_st0_tag == TW_NaN) || (st1_tag == TW_NaN) )
 685     {
 686       real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr);
 687       pop();
 688       return;
 689     }
 690   else if ( (FPU_st0_tag <= TW_Zero) && (st1_tag <= TW_Zero) )
 691     {
 692       /* one of the args is zero, the other valid, or both zero */
 693       if ( FPU_st0_tag == TW_Zero )
 694         {
 695           pop(); FPU_st0_ptr = &st(0);
 696           if ( FPU_st0_ptr->tag == TW_Zero )
 697             arith_invalid(FPU_st0_ptr);   /* Both args zero is invalid */
 698 #ifdef PECULIAR_486
 699           /* This case is not specifically covered in the manual,
 700              but divide-by-zero would seem to be the best response.
 701              However, a real 80486 does it this way... */
 702           else if ( FPU_st0_ptr->tag == TW_Infinity )
 703             {
 704               reg_move(&CONST_INF, FPU_st0_ptr);
 705               return;
 706             }
 707 #endif PECULIAR_486
 708           else
 709             divide_by_zero(st1_ptr->sign^SIGN_NEG^SIGN_POS, FPU_st0_ptr);
 710           return;
 711         }
 712       else
 713         {
 714           /* st(1) contains zero, st(0) valid <> 0 */
 715           /* Zero is the valid answer */
 716           char sign = st1_ptr->sign;
 717 
 718 #ifdef DENORM_OPERAND
 719           if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
 720             return;
 721 #endif DENORM_OPERAND
 722           if ( FPU_st0_ptr->sign == SIGN_NEG )
 723             {
 724               pop(); FPU_st0_ptr = &st(0);
 725               arith_invalid(FPU_st0_ptr);   /* log(negative) */
 726               return;
 727             }
 728 
 729           if ( FPU_st0_ptr->exp < EXP_BIAS ) sign ^= SIGN_NEG^SIGN_POS;
 730           pop(); FPU_st0_ptr = &st(0);
 731           reg_move(&CONST_Z, FPU_st0_ptr);
 732           FPU_st0_ptr->sign = sign;
 733           return;
 734         }
 735     }
 736   /* One or both arg must be an infinity */
 737   else if ( FPU_st0_tag == TW_Infinity )
 738     {
 739       if ( (FPU_st0_ptr->sign == SIGN_NEG) || (st1_tag == TW_Zero) )
 740         { pop(); FPU_st0_ptr = &st(0);
 741           arith_invalid(FPU_st0_ptr);  /* log(-infinity) or 0*log(infinity) */
 742           return; }
 743       else
 744         {
 745           char sign = st1_ptr->sign;
 746 
 747 #ifdef DENORM_OPERAND
 748           if ( (st1_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
 749             return;
 750 #endif DENORM_OPERAND
 751 
 752           pop(); FPU_st0_ptr = &st(0);
 753           reg_move(&CONST_INF, FPU_st0_ptr);
 754           FPU_st0_ptr->sign = sign;
 755           return;
 756         }
 757     }
 758   /* st(1) must be infinity here */
 759   else if ( (FPU_st0_tag == TW_Valid) && (FPU_st0_ptr->sign == SIGN_POS) )
 760     {
 761       if ( FPU_st0_ptr->exp >= EXP_BIAS )
 762         {
 763           if ( (FPU_st0_ptr->exp == EXP_BIAS) &&
 764               (FPU_st0_ptr->sigh == 0x80000000) &&
 765               (FPU_st0_ptr->sigl == 0) )
 766             {
 767               /* st(0) holds 1.0 */
 768               pop(); FPU_st0_ptr = &st(0);
 769               arith_invalid(FPU_st0_ptr); /* infinity*log(1) */
 770               return;
 771             }
 772           /* st(0) is positive and > 1.0 */
 773           pop();
 774         }
 775       else
 776         {
 777           /* st(0) is positive and < 1.0 */
 778 
 779 #ifdef DENORM_OPERAND
 780           if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
 781             return;
 782 #endif DENORM_OPERAND
 783 
 784           st1_ptr->sign ^= SIGN_NEG;
 785           pop();
 786         }
 787       return;
 788     }
 789   else
 790     {
 791       /* st(0) must be zero or negative */
 792       if ( FPU_st0_ptr->tag == TW_Zero )
 793         {
 794           pop(); FPU_st0_ptr = st1_ptr;
 795           st1_ptr->sign ^= SIGN_NEG^SIGN_POS;
 796           /* This should be invalid, but a real 80486 is happy with it. */
 797 #ifndef PECULIAR_486
 798           divide_by_zero(st1_ptr->sign, FPU_st0_ptr);
 799 #endif PECULIAR_486
 800         }
 801       else
 802         {
 803           pop(); FPU_st0_ptr = st1_ptr;
 804           arith_invalid(FPU_st0_ptr);    /* log(negative) */
 805         }
 806       return;
 807     }
 808 }
 809 
 810 
 811 static void fpatan(void)
     /*  */
 812 {
 813   FPU_REG *st1_ptr = &st(1);
 814   char st1_tag = st1_ptr->tag;
 815 
 816   if ( !((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid)) )
 817     {
 818       int saved_control, saved_status;
 819       FPU_REG sum;
 820       int quadrant = st1_ptr->sign | ((FPU_st0_ptr->sign)<<1);
 821 
 822 #ifdef DENORM_OPERAND
 823       if ( ((FPU_st0_ptr->exp <= EXP_UNDER) ||
 824             (st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()) )
 825         return;
 826 #endif DENORM_OPERAND
 827 
 828       /* We use the general purpose arithmetic so we need to save these. */
 829       saved_status = status_word;
 830       saved_control = control_word;
 831       control_word = FULL_PRECISION;
 832 
 833       st1_ptr->sign = FPU_st0_ptr->sign = SIGN_POS;
 834       if ( compare(st1_ptr) == COMP_A_lt_B )
 835         {
 836           quadrant |= 4;
 837           reg_div(FPU_st0_ptr, st1_ptr, &sum, FULL_PRECISION);
 838         }
 839       else
 840         reg_div(st1_ptr, FPU_st0_ptr, &sum, FULL_PRECISION);
 841 
 842       poly_atan(&sum);
 843       
 844       if (quadrant & 4)
 845         {
 846           reg_sub(&CONST_PI2, &sum, &sum, FULL_PRECISION);
 847         }
 848       if (quadrant & 2)
 849         {
 850           reg_sub(&CONST_PI, &sum, &sum, FULL_PRECISION);
 851         }
 852       if (quadrant & 1)
 853         sum.sign ^= SIGN_POS^SIGN_NEG;
 854 
 855       /* All of the basic arithmetic is done now */
 856       control_word = saved_control;
 857       status_word = saved_status;
 858 
 859       reg_move(&sum, st1_ptr);
 860     }
 861   else if ( (FPU_st0_tag == TW_Empty) || (st1_tag == TW_Empty) )
 862     {
 863       stack_underflow_pop(1);
 864       return;
 865     }
 866   else if ( (FPU_st0_tag == TW_NaN) || (st1_tag == TW_NaN) )
 867     {
 868       real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr);
 869       pop();
 870       return;
 871     }
 872   else if ( (FPU_st0_tag == TW_Infinity) || (st1_tag == TW_Infinity) )
 873     {
 874       char sign = st1_ptr->sign;
 875       if ( FPU_st0_tag == TW_Infinity )
 876         {
 877           if ( st1_tag == TW_Infinity )
 878             {
 879               if ( FPU_st0_ptr->sign == SIGN_POS )
 880                 { reg_move(&CONST_PI4, st1_ptr); }
 881               else
 882                 reg_add(&CONST_PI4, &CONST_PI2, st1_ptr, FULL_PRECISION);
 883             }
 884           else
 885             {
 886 
 887 #ifdef DENORM_OPERAND
 888               if ( (st1_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
 889                 return;
 890 #endif DENORM_OPERAND
 891 
 892               if ( FPU_st0_ptr->sign == SIGN_POS )
 893                 { reg_move(&CONST_Z, st1_ptr); pop(); return; }
 894               else
 895                 reg_move(&CONST_PI, st1_ptr);
 896             }
 897         }
 898       else
 899         {
 900           /* st(1) is infinity, st(0) not infinity */
 901 #ifdef DENORM_OPERAND
 902           if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
 903             return;
 904 #endif DENORM_OPERAND
 905 
 906           reg_move(&CONST_PI2, st1_ptr);
 907         }
 908       st1_ptr->sign = sign;
 909     }
 910   else if ( st1_tag == TW_Zero )
 911     {
 912       /* st(0) must be valid or zero */
 913       char sign = st1_ptr->sign;
 914 
 915 #ifdef DENORM_OPERAND
 916       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
 917         return;
 918 #endif DENORM_OPERAND
 919 
 920       if ( FPU_st0_ptr->sign == SIGN_POS )
 921         { reg_move(&CONST_Z, st1_ptr); pop(); return; }
 922       else
 923         reg_move(&CONST_PI, st1_ptr);
 924       st1_ptr->sign = sign;
 925     }
 926   else if ( FPU_st0_tag == TW_Zero )
 927     {
 928       /* st(1) must be TW_Valid here */
 929       char sign = st1_ptr->sign;
 930 
 931 #ifdef DENORM_OPERAND
 932       if ( (st1_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
 933         return;
 934 #endif DENORM_OPERAND
 935 
 936       reg_move(&CONST_PI2, st1_ptr);
 937       st1_ptr->sign = sign;
 938     }
 939 #ifdef PARANOID
 940   else
 941     EXCEPTION(EX_INTERNAL | 0x220);
 942 #endif PARANOID
 943 
 944   pop();
 945   set_precision_flag_up();  /* We do not really know if up or down */
 946 }
 947 
 948 
 949 static void fprem(void)
     /*  */
 950 {
 951   fprem_kernel(RC_CHOP);
 952 }
 953 
 954 
 955 static void fprem1(void)
     /*  */
 956 {
 957   fprem_kernel(RC_RND);
 958 }
 959 
 960 
 961 static void fyl2xp1(void)
     /*  */
 962 {
 963   FPU_REG *st1_ptr = &st(1);
 964   char st1_tag = st1_ptr->tag;
 965 
 966   if ( !((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid)) )
 967     {
 968       int saved_control, saved_status;
 969 
 970 #ifdef DENORM_OPERAND
 971       if ( ((FPU_st0_ptr->exp <= EXP_UNDER) ||
 972             (st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()) )
 973         return;
 974 #endif DENORM_OPERAND
 975 
 976       /* We use the general purpose arithmetic so we need to save these. */
 977       saved_status = status_word;
 978       saved_control = control_word;
 979       control_word = FULL_PRECISION;
 980 
 981       if ( poly_l2p1(FPU_st0_ptr, FPU_st0_ptr) )
 982         {
 983           arith_invalid(st1_ptr);  /* poly_l2p1() returned invalid */
 984           pop(); return;
 985         }
 986       
 987       /* Enough of the basic arithmetic is done now */
 988       control_word = saved_control;
 989       status_word = saved_status;
 990 
 991       /* Let the multiply set the flags */
 992       reg_mul(FPU_st0_ptr, st1_ptr, st1_ptr, FULL_PRECISION);
 993 
 994       pop();
 995     }
 996   else if ( (FPU_st0_tag == TW_Empty) | (st1_tag == TW_Empty) )
 997     {
 998       stack_underflow_pop(1);
 999       return;
1000     }
1001   else if ( FPU_st0_tag == TW_Zero )
1002     {
1003       if ( st1_tag <= TW_Zero )
1004         {
1005 
1006 #ifdef DENORM_OPERAND
1007               if ( (st1_tag == TW_Valid) && (st1_ptr->exp <= EXP_UNDER) &&
1008                   (denormal_operand()) )
1009                 return;
1010 #endif DENORM_OPERAND
1011 
1012           st1_ptr->sign ^= FPU_st0_ptr->sign;
1013           reg_move(FPU_st0_ptr, st1_ptr);
1014         }
1015       else if ( st1_tag == TW_Infinity )
1016         {
1017           arith_invalid(st1_ptr);  /* Infinity*log(1) */
1018           pop();
1019           return;
1020         }
1021       else if ( st1_tag == TW_NaN )
1022         {
1023           real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr);
1024           pop();
1025           return;
1026         }
1027 #ifdef PARANOID
1028       else
1029         {
1030           EXCEPTION(EX_INTERNAL | 0x116);
1031           return;
1032         }
1033 #endif PARANOID
1034       pop(); return;
1035     }
1036   else if ( FPU_st0_tag == TW_Valid )
1037     {
1038       if ( st1_tag == TW_Zero )
1039         {
1040           if ( FPU_st0_ptr->sign == SIGN_NEG )
1041             {
1042               if ( FPU_st0_ptr->exp >= EXP_BIAS )
1043                 {
1044                   /* st(0) holds <= -1.0 */
1045                   arith_invalid(st1_ptr); /* infinity*log(1) */
1046                   pop(); return;
1047                 }
1048 #ifdef DENORM_OPERAND
1049               if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
1050                 return;
1051 #endif DENORM_OPERAND
1052               st1_ptr->sign ^= SIGN_POS^SIGN_NEG;
1053               pop(); return;
1054             }
1055 #ifdef DENORM_OPERAND
1056           if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
1057             return;
1058 #endif DENORM_OPERAND
1059           pop(); return;
1060         }
1061       if ( st1_tag == TW_Infinity )
1062         {
1063           if ( FPU_st0_ptr->sign == SIGN_NEG )
1064             {
1065               if ( (FPU_st0_ptr->exp >= EXP_BIAS) &&
1066                   !((FPU_st0_ptr->sigh == 0x80000000) &&
1067                     (FPU_st0_ptr->sigl == 0)) )
1068                 {
1069                   /* st(0) holds < -1.0 */
1070                   arith_invalid(st1_ptr);
1071                   pop(); return;
1072                 }
1073 #ifdef DENORM_OPERAND
1074               if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
1075                 return;
1076 #endif DENORM_OPERAND
1077               st1_ptr->sign ^= SIGN_POS^SIGN_NEG;
1078               pop(); return;
1079             }
1080 #ifdef DENORM_OPERAND
1081           if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
1082             return;
1083 #endif DENORM_OPERAND
1084           pop(); return;
1085         }
1086       if ( st1_tag == TW_NaN )
1087         {
1088           real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr);
1089           pop();
1090           return;
1091         }
1092     }
1093   else if ( FPU_st0_tag == TW_NaN )
1094     {
1095       real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr);
1096       pop();
1097       return;
1098     }
1099   else if ( FPU_st0_tag == TW_Infinity )
1100     {
1101       if ( st1_tag == TW_NaN )
1102         {
1103           real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr);
1104           pop();
1105           return;
1106         }
1107       else if ( (FPU_st0_ptr->sign == SIGN_NEG) ||
1108                (st1_tag == TW_Zero) )
1109         {
1110           arith_invalid(st1_ptr);  /* log(infinity) */
1111           pop();
1112           return;
1113         }
1114         
1115       /* st(1) must be valid here. */
1116 
1117 #ifdef DENORM_OPERAND
1118           if ( (st1_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
1119             return;
1120 #endif DENORM_OPERAND
1121 
1122       /* The Manual says that log(Infinity) is invalid, but a real
1123          80486 sensibly says that it is o.k. */
1124       { char sign = st1_ptr->sign;
1125         reg_move(&CONST_INF, st1_ptr);
1126         st1_ptr->sign = sign;
1127       }
1128       pop();
1129       return;
1130     }
1131 #ifdef PARANOID
1132   else
1133     {
1134       EXCEPTION(EX_INTERNAL | 0x117);
1135     }
1136 #endif PARANOID
1137 }
1138 
1139 
1140 static void fscale(void)
     /*  */
1141 {
1142   FPU_REG *st1_ptr = &st(1);
1143   char st1_tag = st1_ptr->tag;
1144   int old_cw = control_word;
1145 
1146   if ( !((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid)) )
1147     {
1148       long scale;
1149       FPU_REG tmp;
1150 
1151 #ifdef DENORM_OPERAND
1152       if ( ((FPU_st0_ptr->exp <= EXP_UNDER) ||
1153             (st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()) )
1154         return;
1155 #endif DENORM_OPERAND
1156 
1157       if ( st1_ptr->exp > EXP_BIAS + 30 )
1158         {
1159           /* 2^31 is far too large, would require 2^(2^30) or 2^(-2^30) */
1160           char sign;
1161 
1162           if ( st1_ptr->sign == SIGN_POS )
1163             {
1164               EXCEPTION(EX_Overflow);
1165               sign = FPU_st0_ptr->sign;
1166               reg_move(&CONST_INF, FPU_st0_ptr);
1167               FPU_st0_ptr->sign = sign;
1168             }
1169           else
1170             {
1171               EXCEPTION(EX_Underflow);
1172               sign = FPU_st0_ptr->sign;
1173               reg_move(&CONST_Z, FPU_st0_ptr);
1174               FPU_st0_ptr->sign = sign;
1175             }
1176           return;
1177         }
1178 
1179       control_word &= ~CW_RC;
1180       control_word |= RC_CHOP;
1181       reg_move(st1_ptr, &tmp);
1182       round_to_int(&tmp);               /* This can never overflow here */
1183       control_word = old_cw;
1184       scale = st1_ptr->sign ? -tmp.sigl : tmp.sigl;
1185       scale += FPU_st0_ptr->exp;
1186       FPU_st0_ptr->exp = scale;
1187 
1188       /* Use round_reg() to properly detect under/overflow etc */
1189       round_reg(FPU_st0_ptr, 0, control_word);
1190 
1191       return;
1192     }
1193   else if ( FPU_st0_tag == TW_Valid )
1194     {
1195       if ( st1_tag == TW_Zero )
1196         {
1197 
1198 #ifdef DENORM_OPERAND
1199           if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
1200             return;
1201 #endif DENORM_OPERAND
1202 
1203           return;
1204         }
1205       if ( st1_tag == TW_Infinity )
1206         {
1207           char sign = st1_ptr->sign;
1208 
1209 #ifdef DENORM_OPERAND
1210           if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
1211             return;
1212 #endif DENORM_OPERAND
1213 
1214           if ( sign == SIGN_POS )
1215             { reg_move(&CONST_INF, FPU_st0_ptr); }
1216           else
1217               reg_move(&CONST_Z, FPU_st0_ptr);
1218           FPU_st0_ptr->sign = sign;
1219           return;
1220         }
1221       if ( st1_tag == TW_NaN )
1222         { real_2op_NaN(FPU_st0_ptr, st1_ptr, FPU_st0_ptr); return; }
1223     }
1224   else if ( FPU_st0_tag == TW_Zero )
1225     {
1226       if ( st1_tag == TW_Valid )
1227         {
1228 
1229 #ifdef DENORM_OPERAND
1230           if ( (st1_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
1231             return;
1232 #endif DENORM_OPERAND
1233 
1234           return;
1235         }
1236       else if ( st1_tag == TW_Zero ) { return; }
1237       else if ( st1_tag == TW_Infinity )
1238         {
1239           if ( st1_ptr->sign == SIGN_NEG )
1240             return;
1241           else
1242             {
1243               arith_invalid(FPU_st0_ptr); /* Zero scaled by +Infinity */
1244               return;
1245             }
1246         }
1247       else if ( st1_tag == TW_NaN )
1248         { real_2op_NaN(FPU_st0_ptr, st1_ptr, FPU_st0_ptr); return; }
1249     }
1250   else if ( FPU_st0_tag == TW_Infinity )
1251     {
1252       if ( st1_tag == TW_Valid )
1253         {
1254 
1255 #ifdef DENORM_OPERAND
1256           if ( (st1_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
1257             return;
1258 #endif DENORM_OPERAND
1259 
1260           return;
1261         }
1262       if ( ((st1_tag == TW_Infinity) && (st1_ptr->sign == SIGN_POS))
1263           || (st1_tag == TW_Zero) )
1264         return;
1265       else if ( st1_tag == TW_Infinity )
1266         {
1267           arith_invalid(FPU_st0_ptr); /* Infinity scaled by -Infinity */
1268           return;
1269         }
1270       else if ( st1_tag == TW_NaN )
1271         { real_2op_NaN(FPU_st0_ptr, st1_ptr, FPU_st0_ptr); return; }
1272     }
1273   else if ( FPU_st0_tag == TW_NaN )
1274     {
1275       if ( st1_tag != TW_Empty )
1276         { real_2op_NaN(FPU_st0_ptr, st1_ptr, FPU_st0_ptr); return; }
1277     }
1278 
1279 #ifdef PARANOID
1280   if ( !((FPU_st0_tag == TW_Empty) || (st1_tag == TW_Empty)) )
1281     {
1282       EXCEPTION(EX_INTERNAL | 0x115);
1283       return;
1284     }
1285 #endif
1286 
1287   /* At least one of st(0), st(1) must be empty */
1288   stack_underflow();
1289 
1290 }
1291 
1292 
1293 /*---------------------------------------------------------------------------*/
1294 
1295 static FUNC trig_table_a[] = {
1296   f2xm1, fyl2x, fptan, fpatan, fxtract, fprem1, fdecstp, fincstp
1297 };
1298 
1299 void trig_a(void)
     /*  */
1300 {
1301   (trig_table_a[FPU_rm])();
1302 }
1303 
1304 
1305 static FUNC trig_table_b[] =
1306   {
1307     fprem, fyl2xp1, fsqrt_, fsincos, frndint_, fscale, fsin, fcos
1308   };
1309 
1310 void trig_b(void)
     /*  */
1311 {
1312   (trig_table_b[FPU_rm])();
1313 }
/* */
root/kernel/FPU-emu/fpu_trig.c

DEFINITIONS
