PageRenderTime 129ms CodeModel.GetById 17ms app.highlight 101ms RepoModel.GetById 1ms app.codeStats 1ms

/arch/i386/math-emu/fpu_trig.c

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