arch/x86/math-emu/fpu_trig.c

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