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 fsin(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 int f_cos(FPU_REG *st0_ptr, u_char tag)
 612 {
 613         u_char st0_sign;
 614
 615         st0_sign = getsign(st0_ptr);
 616
 617         if (tag == TAG_Valid) {
 618                 int q;
 619
 620                 if (exponent(st0_ptr) > -40) {
 621                         if ((exponent(st0_ptr) < 0)
 622                             || ((exponent(st0_ptr) == 0)
 623                                 && (significand(st0_ptr) <=
 624                                     0xc90fdaa22168c234LL))) {
 625                                 poly_cos(st0_ptr);
 626
 627                                 /* We do not really know if up or down */
 628                                 set_precision_flag_down();
 629
 630                                 return 0;
 631                         } else if ((q = trig_arg(st0_ptr, FCOS)) != -1) {
 632                                 poly_sine(st0_ptr);
 633
 634                                 if ((q + 1) & 2)
 635                                         changesign(st0_ptr);
 636
 637                                 /* We do not really know if up or down */
 638                                 set_precision_flag_down();
 639
 640                                 return 0;
 641                         } else {
 642                                 /* Operand is out of range */
 643                                 return 1;
 644                         }
 645                 } else {
 646                       denormal_arg:
 647
 648                         setcc(0);
 649                         FPU_copy_to_reg0(&CONST_1, TAG_Valid);
 650 #ifdef PECULIAR_486
 651                         set_precision_flag_down();      /* 80486 appears to do this. */
 652 #else
 653                         set_precision_flag_up();        /* Must be up. */
 654 #endif /* PECULIAR_486 */
 655                         return 0;
 656                 }
 657         } else if (tag == TAG_Zero) {
 658                 FPU_copy_to_reg0(&CONST_1, TAG_Valid);
 659                 setcc(0);
 660                 return 0;
 661         }
 662
 663         if (tag == TAG_Special)
 664                 tag = FPU_Special(st0_ptr);
 665
 666         if (tag == TW_Denormal) {
 667                 if (denormal_operand() < 0)
 668                         return 1;
 669
 670                 goto denormal_arg;
 671         } else if (tag == TW_Infinity) {
 672                 /* The 80486 treats infinity as an invalid operand */
 673                 arith_invalid(0);
 674                 return 1;
 675         } else {
 676                 single_arg_error(st0_ptr, tag); /* requires st0_ptr == &st(0) */
 677                 return 1;
 678         }
 679 }
 680
 681 static void fcos(FPU_REG *st0_ptr, u_char st0_tag)
 682 {
 683         f_cos(st0_ptr, st0_tag);
 684 }
 685
 686 static void fsincos(FPU_REG *st0_ptr, u_char st0_tag)
 687 {
 688         FPU_REG *st_new_ptr;
 689         FPU_REG arg;
 690         u_char tag;
 691
 692         /* Stack underflow has higher priority */
 693         if (st0_tag == TAG_Empty) {
 694                 FPU_stack_underflow();  /* Puts a QNaN in st(0) */
 695                 if (control_word & CW_Invalid) {
 696                         st_new_ptr = &st(-1);
 697                         push();
 698                         FPU_stack_underflow();  /* Puts a QNaN in the new st(0) */
 699                 }
 700                 return;
 701         }
 702
 703         if (STACK_OVERFLOW) {
 704                 FPU_stack_overflow();
 705                 return;
 706         }
 707
 708         if (st0_tag == TAG_Special)
 709                 tag = FPU_Special(st0_ptr);
 710         else
 711                 tag = st0_tag;
 712
 713         if (tag == TW_NaN) {
 714                 single_arg_2_error(st0_ptr, TW_NaN);
 715                 return;
 716         } else if (tag == TW_Infinity) {
 717                 /* The 80486 treats infinity as an invalid operand */
 718                 if (arith_invalid(0) >= 0) {
 719                         /* Masked response */
 720                         push();
 721                         arith_invalid(0);
 722                 }
 723                 return;
 724         }
 725
 726         reg_copy(st0_ptr, &arg);
 727         if (!fsin(st0_ptr, st0_tag)) {
 728                 push();
 729                 FPU_copy_to_reg0(&arg, st0_tag);
 730                 f_cos(&st(0), st0_tag);
 731         } else {
 732                 /* An error, so restore st(0) */
 733                 FPU_copy_to_reg0(&arg, st0_tag);
 734         }
 735 }
 736
 737 /*---------------------------------------------------------------------------*/
 738 /* The following all require two arguments: st(0) and st(1) */
 739
 740 /* A lean, mean kernel for the fprem instructions. This relies upon
 741    the division and rounding to an integer in do_fprem giving an
 742    exact result. Because of this, rem_kernel() needs to deal only with
 743    the least significant 64 bits, the more significant bits of the
 744    result must be zero.
 745  */
 746 static void rem_kernel(unsigned long long st0, unsigned long long *y,
 747                        unsigned long long st1, unsigned long long q, int n)
 748 {
 749         int dummy;
 750         unsigned long long x;
 751
 752         x = st0 << n;
 753
 754         /* Do the required multiplication and subtraction in the one operation */
 755
 756         /* lsw x -= lsw st1 * lsw q */
 757         asm volatile ("mull %4; subl %%eax,%0; sbbl %%edx,%1":"=m"
 758                       (((unsigned *)&x)[0]), "=m"(((unsigned *)&x)[1]),
 759                       "=a"(dummy)
 760                       :"2"(((unsigned *)&st1)[0]), "m"(((unsigned *)&q)[0])
 761                       :"%dx");
 762         /* msw x -= msw st1 * lsw q */
 763         asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x)[1]),
 764                       "=a"(dummy)
 765                       :"1"(((unsigned *)&st1)[1]), "m"(((unsigned *)&q)[0])
 766                       :"%dx");
 767         /* msw x -= lsw st1 * msw q */
 768         asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x)[1]),
 769                       "=a"(dummy)
 770                       :"1"(((unsigned *)&st1)[0]), "m"(((unsigned *)&q)[1])
 771                       :"%dx");
 772
 773         *y = x;
 774 }
 775
 776 /* Remainder of st(0) / st(1) */
 777 /* This routine produces exact results, i.e. there is never any
 778    rounding or truncation, etc of the result. */
 779 static void do_fprem(FPU_REG *st0_ptr, u_char st0_tag, int round)
 780 {
 781         FPU_REG *st1_ptr = &st(1);
 782         u_char st1_tag = FPU_gettagi(1);
 783
 784         if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
 785                 FPU_REG tmp, st0, st1;
 786                 u_char st0_sign, st1_sign;
 787                 u_char tmptag;
 788                 int tag;
 789                 int old_cw;
 790                 int expdif;
 791                 long long q;
 792                 unsigned short saved_status;
 793                 int cc;
 794
 795               fprem_valid:
 796                 /* Convert registers for internal use. */
 797                 st0_sign = FPU_to_exp16(st0_ptr, &st0);
 798                 st1_sign = FPU_to_exp16(st1_ptr, &st1);
 799                 expdif = exponent16(&st0) - exponent16(&st1);
 800
 801                 old_cw = control_word;
 802                 cc = 0;
 803
 804                 /* We want the status following the denorm tests, but don't want
 805                    the status changed by the arithmetic operations. */
 806                 saved_status = partial_status;
 807                 control_word &= ~CW_RC;
 808                 control_word |= RC_CHOP;
 809
 810                 if (expdif < 64) {
 811                         /* This should be the most common case */
 812
 813                         if (expdif > -2) {
 814                                 u_char sign = st0_sign ^ st1_sign;
 815                                 tag = FPU_u_div(&st0, &st1, &tmp,
 816                                                 PR_64_BITS | RC_CHOP | 0x3f,
 817                                                 sign);
 818                                 setsign(&tmp, sign);
 819
 820                                 if (exponent(&tmp) >= 0) {
 821                                         FPU_round_to_int(&tmp, tag);    /* Fortunately, this can't
 822                                                                            overflow to 2^64 */
 823                                         q = significand(&tmp);
 824
 825                                         rem_kernel(significand(&st0),
 826                                                    &significand(&tmp),
 827                                                    significand(&st1),
 828                                                    q, expdif);
 829
 830                                         setexponent16(&tmp, exponent16(&st1));
 831                                 } else {
 832                                         reg_copy(&st0, &tmp);
 833                                         q = 0;
 834                                 }
 835
 836                                 if ((round == RC_RND)
 837                                     && (tmp.sigh & 0xc0000000)) {
 838                                         /* We may need to subtract st(1) once more,
 839                                            to get a result <= 1/2 of st(1). */
 840                                         unsigned long long x;
 841                                         expdif =
 842                                             exponent16(&st1) - exponent16(&tmp);
 843                                         if (expdif <= 1) {
 844                                                 if (expdif == 0)
 845                                                         x = significand(&st1) -
 846                                                             significand(&tmp);
 847                                                 else    /* expdif is 1 */
 848                                                         x = (significand(&st1)
 849                                                              << 1) -
 850                                                             significand(&tmp);
 851                                                 if ((x < significand(&tmp)) ||
 852                                                     /* or equi-distant (from 0 & st(1)) and q is odd */
 853                                                     ((x == significand(&tmp))
 854                                                      && (q & 1))) {
 855                                                         st0_sign = !st0_sign;
 856                                                         significand(&tmp) = x;
 857                                                         q++;
 858                                                 }
 859                                         }
 860                                 }
 861
 862                                 if (q & 4)
 863                                         cc |= SW_C0;
 864                                 if (q & 2)
 865                                         cc |= SW_C3;
 866                                 if (q & 1)
 867                                         cc |= SW_C1;
 868                         } else {
 869                                 control_word = old_cw;
 870                                 setcc(0);
 871                                 return;
 872                         }
 873                 } else {
 874                         /* There is a large exponent difference ( >= 64 ) */
 875                         /* To make much sense, the code in this section should
 876                            be done at high precision. */
 877                         int exp_1, N;
 878                         u_char sign;
 879
 880                         /* prevent overflow here */
 881                         /* N is 'a number between 32 and 63' (p26-113) */
 882                         reg_copy(&st0, &tmp);
 883                         tmptag = st0_tag;
 884                         N = (expdif & 0x0000001f) + 32; /* This choice gives results
 885                                                            identical to an AMD 486 */
 886                         setexponent16(&tmp, N);
 887                         exp_1 = exponent16(&st1);
 888                         setexponent16(&st1, 0);
 889                         expdif -= N;
 890
 891                         sign = getsign(&tmp) ^ st1_sign;
 892                         tag =
 893                             FPU_u_div(&tmp, &st1, &tmp,
 894                                       PR_64_BITS | RC_CHOP | 0x3f, sign);
 895                         setsign(&tmp, sign);
 896
 897                         FPU_round_to_int(&tmp, tag);    /* Fortunately, this can't
 898                                                            overflow to 2^64 */
 899
 900                         rem_kernel(significand(&st0),
 901                                    &significand(&tmp),
 902                                    significand(&st1),
 903                                    significand(&tmp), exponent(&tmp)
 904                             );
 905                         setexponent16(&tmp, exp_1 + expdif);
 906
 907                         /* It is possible for the operation to be complete here.
 908                            What does the IEEE standard say? The Intel 80486 manual
 909                            implies that the operation will never be completed at this
 910                            point, and the behaviour of a real 80486 confirms this.
 911                          */
 912                         if (!(tmp.sigh | tmp.sigl)) {
 913                                 /* The result is zero */
 914                                 control_word = old_cw;
 915                                 partial_status = saved_status;
 916                                 FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
 917                                 setsign(&st0, st0_sign);
 918 #ifdef PECULIAR_486
 919                                 setcc(SW_C2);
 920 #else
 921                                 setcc(0);
 922 #endif /* PECULIAR_486 */
 923                                 return;
 924                         }
 925                         cc = SW_C2;
 926                 }
 927
 928                 control_word = old_cw;
 929                 partial_status = saved_status;
 930                 tag = FPU_normalize_nuo(&tmp);
 931                 reg_copy(&tmp, st0_ptr);
 932
 933                 /* The only condition to be looked for is underflow,
 934                    and it can occur here only if underflow is unmasked. */
 935                 if ((exponent16(&tmp) <= EXP_UNDER) && (tag != TAG_Zero)
 936                     && !(control_word & CW_Underflow)) {
 937                         setcc(cc);
 938                         tag = arith_underflow(st0_ptr);
 939                         setsign(st0_ptr, st0_sign);
 940                         FPU_settag0(tag);
 941                         return;
 942                 } else if ((exponent16(&tmp) > EXP_UNDER) || (tag == TAG_Zero)) {
 943                         stdexp(st0_ptr);
 944                         setsign(st0_ptr, st0_sign);
 945                 } else {
 946                         tag =
 947                             FPU_round(st0_ptr, 0, 0, FULL_PRECISION, st0_sign);
 948                 }
 949                 FPU_settag0(tag);
 950                 setcc(cc);
 951
 952                 return;
 953         }
 954
 955         if (st0_tag == TAG_Special)
 956                 st0_tag = FPU_Special(st0_ptr);
 957         if (st1_tag == TAG_Special)
 958                 st1_tag = FPU_Special(st1_ptr);
 959
 960         if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
 961             || ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
 962             || ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
 963                 if (denormal_operand() < 0)
 964                         return;
 965                 goto fprem_valid;
 966         } else if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
 967                 FPU_stack_underflow();
 968                 return;
 969         } else if (st0_tag == TAG_Zero) {
 970                 if (st1_tag == TAG_Valid) {
 971                         setcc(0);
 972                         return;
 973                 } else if (st1_tag == TW_Denormal) {
 974                         if (denormal_operand() < 0)
 975                                 return;
 976                         setcc(0);
 977                         return;
 978                 } else if (st1_tag == TAG_Zero) {
 979                         arith_invalid(0);
 980                         return;
 981                 } /* fprem(?,0) always invalid */
 982                 else if (st1_tag == TW_Infinity) {
 983                         setcc(0);
 984                         return;
 985                 }
 986         } else if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
 987                 if (st1_tag == TAG_Zero) {
 988                         arith_invalid(0);       /* fprem(Valid,Zero) is invalid */
 989                         return;
 990                 } else if (st1_tag != TW_NaN) {
 991                         if (((st0_tag == TW_Denormal)
 992                              || (st1_tag == TW_Denormal))
 993                             && (denormal_operand() < 0))
 994                                 return;
 995
 996                         if (st1_tag == TW_Infinity) {
 997                                 /* fprem(Valid,Infinity) is o.k. */
 998                                 setcc(0);
 999                                 return;
1000                         }
1001                 }
1002         } else if (st0_tag == TW_Infinity) {
1003                 if (st1_tag != TW_NaN) {
1004                         arith_invalid(0);       /* fprem(Infinity,?) is invalid */
1005                         return;
1006                 }
1007         }
1008
1009         /* One of the registers must contain a NaN if we got here. */
1010
1011 #ifdef PARANOID
1012         if ((st0_tag != TW_NaN) && (st1_tag != TW_NaN))
1013                 EXCEPTION(EX_INTERNAL | 0x118);
1014 #endif /* PARANOID */
1015
1016         real_2op_NaN(st1_ptr, st1_tag, 0, st1_ptr);
1017
1018 }
1019
1020 /* ST(1) <- ST(1) * log ST;  pop ST */
1021 static void fyl2x(FPU_REG *st0_ptr, u_char st0_tag)
1022 {
1023         FPU_REG *st1_ptr = &st(1), exponent;
1024         u_char st1_tag = FPU_gettagi(1);
1025         u_char sign;
1026         int e, tag;
1027
1028         clear_C1();
1029
1030         if ((st0_tag == TAG_Valid) && (st1_tag == TAG_Valid)) {
1031               both_valid:
1032                 /* Both regs are Valid or Denormal */
1033                 if (signpositive(st0_ptr)) {
1034                         if (st0_tag == TW_Denormal)
1035                                 FPU_to_exp16(st0_ptr, st0_ptr);
1036                         else
1037                                 /* Convert st(0) for internal use. */
1038                                 setexponent16(st0_ptr, exponent(st0_ptr));
1039
1040                         if ((st0_ptr->sigh == 0x80000000)
1041                             && (st0_ptr->sigl == 0)) {
1042                                 /* Special case. The result can be precise. */
1043                                 u_char esign;
1044                                 e = exponent16(st0_ptr);
1045                                 if (e >= 0) {
1046                                         exponent.sigh = e;
1047                                         esign = SIGN_POS;
1048                                 } else {
1049                                         exponent.sigh = -e;
1050                                         esign = SIGN_NEG;
1051                                 }
1052                                 exponent.sigl = 0;
1053                                 setexponent16(&exponent, 31);
1054                                 tag = FPU_normalize_nuo(&exponent);
1055                                 stdexp(&exponent);
1056                                 setsign(&exponent, esign);
1057                                 tag =
1058                                     FPU_mul(&exponent, tag, 1, FULL_PRECISION);
1059                                 if (tag >= 0)
1060                                         FPU_settagi(1, tag);
1061                         } else {
1062                                 /* The usual case */
1063                                 sign = getsign(st1_ptr);
1064                                 if (st1_tag == TW_Denormal)
1065                                         FPU_to_exp16(st1_ptr, st1_ptr);
1066                                 else
1067                                         /* Convert st(1) for internal use. */
1068                                         setexponent16(st1_ptr,
1069                                                       exponent(st1_ptr));
1070                                 poly_l2(st0_ptr, st1_ptr, sign);
1071                         }
1072                 } else {
1073                         /* negative */
1074                         if (arith_invalid(1) < 0)
1075                                 return;
1076                 }
1077
1078                 FPU_pop();
1079
1080                 return;
1081         }
1082
1083         if (st0_tag == TAG_Special)
1084                 st0_tag = FPU_Special(st0_ptr);
1085         if (st1_tag == TAG_Special)
1086                 st1_tag = FPU_Special(st1_ptr);
1087
1088         if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
1089                 FPU_stack_underflow_pop(1);
1090                 return;
1091         } else if ((st0_tag <= TW_Denormal) && (st1_tag <= TW_Denormal)) {
1092                 if (st0_tag == TAG_Zero) {
1093                         if (st1_tag == TAG_Zero) {
1094                                 /* Both args zero is invalid */
1095                                 if (arith_invalid(1) < 0)
1096                                         return;
1097                         } else {
1098                                 u_char sign;
1099                                 sign = getsign(st1_ptr) ^ SIGN_NEG;
1100                                 if (FPU_divide_by_zero(1, sign) < 0)
1101                                         return;
1102
1103                                 setsign(st1_ptr, sign);
1104                         }
1105                 } else if (st1_tag == TAG_Zero) {
1106                         /* st(1) contains zero, st(0) valid <> 0 */
1107                         /* Zero is the valid answer */
1108                         sign = getsign(st1_ptr);
1109
1110                         if (signnegative(st0_ptr)) {
1111                                 /* log(negative) */
1112                                 if (arith_invalid(1) < 0)
1113                                         return;
1114                         } else if ((st0_tag == TW_Denormal)
1115                                    && (denormal_operand() < 0))
1116                                 return;
1117                         else {
1118                                 if (exponent(st0_ptr) < 0)
1119                                         sign ^= SIGN_NEG;
1120
1121                                 FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
1122                                 setsign(st1_ptr, sign);
1123                         }
1124                 } else {
1125                         /* One or both operands are denormals. */
1126                         if (denormal_operand() < 0)
1127                                 return;
1128                         goto both_valid;
1129                 }
1130         } else if ((st0_tag == TW_NaN) || (st1_tag == TW_NaN)) {
1131                 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1132                         return;
1133         }
1134         /* One or both arg must be an infinity */
1135         else if (st0_tag == TW_Infinity) {
1136                 if ((signnegative(st0_ptr)) || (st1_tag == TAG_Zero)) {
1137                         /* log(-infinity) or 0*log(infinity) */
1138                         if (arith_invalid(1) < 0)
1139                                 return;
1140                 } else {
1141                         u_char sign = getsign(st1_ptr);
1142
1143                         if ((st1_tag == TW_Denormal)
1144                             && (denormal_operand() < 0))
1145                                 return;
1146
1147                         FPU_copy_to_reg1(&CONST_INF, TAG_Special);
1148                         setsign(st1_ptr, sign);
1149                 }
1150         }
1151         /* st(1) must be infinity here */
1152         else if (((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal))
1153                  && (signpositive(st0_ptr))) {
1154                 if (exponent(st0_ptr) >= 0) {
1155                         if ((exponent(st0_ptr) == 0) &&
1156                             (st0_ptr->sigh == 0x80000000) &&
1157                             (st0_ptr->sigl == 0)) {
1158                                 /* st(0) holds 1.0 */
1159                                 /* infinity*log(1) */
1160                                 if (arith_invalid(1) < 0)
1161                                         return;
1162                         }
1163                         /* else st(0) is positive and > 1.0 */
1164                 } else {
1165                         /* st(0) is positive and < 1.0 */
1166
1167                         if ((st0_tag == TW_Denormal)
1168                             && (denormal_operand() < 0))
1169                                 return;
1170
1171                         changesign(st1_ptr);
1172                 }
1173         } else {
1174                 /* st(0) must be zero or negative */
1175                 if (st0_tag == TAG_Zero) {
1176                         /* This should be invalid, but a real 80486 is happy with it. */
1177
1178 #ifndef PECULIAR_486
1179                         sign = getsign(st1_ptr);
1180                         if (FPU_divide_by_zero(1, sign) < 0)
1181                                 return;
1182 #endif /* PECULIAR_486 */
1183
1184                         changesign(st1_ptr);
1185                 } else if (arith_invalid(1) < 0)        /* log(negative) */
1186                         return;
1187         }
1188
1189         FPU_pop();
1190 }
1191
1192 static void fpatan(FPU_REG *st0_ptr, u_char st0_tag)
1193 {
1194         FPU_REG *st1_ptr = &st(1);
1195         u_char st1_tag = FPU_gettagi(1);
1196         int tag;
1197
1198         clear_C1();
1199         if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1200               valid_atan:
1201
1202                 poly_atan(st0_ptr, st0_tag, st1_ptr, st1_tag);
1203
1204                 FPU_pop();
1205
1206                 return;
1207         }
1208
1209         if (st0_tag == TAG_Special)
1210                 st0_tag = FPU_Special(st0_ptr);
1211         if (st1_tag == TAG_Special)
1212                 st1_tag = FPU_Special(st1_ptr);
1213
1214         if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
1215             || ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
1216             || ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
1217                 if (denormal_operand() < 0)
1218                         return;
1219
1220                 goto valid_atan;
1221         } else if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
1222                 FPU_stack_underflow_pop(1);
1223                 return;
1224         } else if ((st0_tag == TW_NaN) || (st1_tag == TW_NaN)) {
1225                 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) >= 0)
1226                         FPU_pop();
1227                 return;
1228         } else if ((st0_tag == TW_Infinity) || (st1_tag == TW_Infinity)) {
1229                 u_char sign = getsign(st1_ptr);
1230                 if (st0_tag == TW_Infinity) {
1231                         if (st1_tag == TW_Infinity) {
1232                                 if (signpositive(st0_ptr)) {
1233                                         FPU_copy_to_reg1(&CONST_PI4, TAG_Valid);
1234                                 } else {
1235                                         setpositive(st1_ptr);
1236                                         tag =
1237                                             FPU_u_add(&CONST_PI4, &CONST_PI2,
1238                                                       st1_ptr, FULL_PRECISION,
1239                                                       SIGN_POS,
1240                                                       exponent(&CONST_PI4),
1241                                                       exponent(&CONST_PI2));
1242                                         if (tag >= 0)
1243                                                 FPU_settagi(1, tag);
1244                                 }
1245                         } else {
1246                                 if ((st1_tag == TW_Denormal)
1247                                     && (denormal_operand() < 0))
1248                                         return;
1249
1250                                 if (signpositive(st0_ptr)) {
1251                                         FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
1252                                         setsign(st1_ptr, sign); /* An 80486 preserves the sign */
1253                                         FPU_pop();
1254                                         return;
1255                                 } else {
1256                                         FPU_copy_to_reg1(&CONST_PI, TAG_Valid);
1257                                 }
1258                         }
1259                 } else {
1260                         /* st(1) is infinity, st(0) not infinity */
1261                         if ((st0_tag == TW_Denormal)
1262                             && (denormal_operand() < 0))
1263                                 return;
1264
1265                         FPU_copy_to_reg1(&CONST_PI2, TAG_Valid);
1266                 }
1267                 setsign(st1_ptr, sign);
1268         } else if (st1_tag == TAG_Zero) {
1269                 /* st(0) must be valid or zero */
1270                 u_char sign = getsign(st1_ptr);
1271
1272                 if ((st0_tag == TW_Denormal) && (denormal_operand() < 0))
1273                         return;
1274
1275                 if (signpositive(st0_ptr)) {
1276                         /* An 80486 preserves the sign */
1277                         FPU_pop();
1278                         return;
1279                 }
1280
1281                 FPU_copy_to_reg1(&CONST_PI, TAG_Valid);
1282                 setsign(st1_ptr, sign);
1283         } else if (st0_tag == TAG_Zero) {
1284                 /* st(1) must be TAG_Valid here */
1285                 u_char sign = getsign(st1_ptr);
1286
1287                 if ((st1_tag == TW_Denormal) && (denormal_operand() < 0))
1288                         return;
1289
1290                 FPU_copy_to_reg1(&CONST_PI2, TAG_Valid);
1291                 setsign(st1_ptr, sign);
1292         }
1293 #ifdef PARANOID
1294         else
1295                 EXCEPTION(EX_INTERNAL | 0x125);
1296 #endif /* PARANOID */
1297
1298         FPU_pop();
1299         set_precision_flag_up();        /* We do not really know if up or down */
1300 }
1301
1302 static void fprem(FPU_REG *st0_ptr, u_char st0_tag)
1303 {
1304         do_fprem(st0_ptr, st0_tag, RC_CHOP);
1305 }
1306
1307 static void fprem1(FPU_REG *st0_ptr, u_char st0_tag)
1308 {
1309         do_fprem(st0_ptr, st0_tag, RC_RND);
1310 }
1311
1312 static void fyl2xp1(FPU_REG *st0_ptr, u_char st0_tag)
1313 {
1314         u_char sign, sign1;
1315         FPU_REG *st1_ptr = &st(1), a, b;
1316         u_char st1_tag = FPU_gettagi(1);
1317
1318         clear_C1();
1319         if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1320               valid_yl2xp1:
1321
1322                 sign = getsign(st0_ptr);
1323                 sign1 = getsign(st1_ptr);
1324
1325                 FPU_to_exp16(st0_ptr, &a);
1326                 FPU_to_exp16(st1_ptr, &b);
1327
1328                 if (poly_l2p1(sign, sign1, &a, &b, st1_ptr))
1329                         return;
1330
1331                 FPU_pop();
1332                 return;
1333         }
1334
1335         if (st0_tag == TAG_Special)
1336                 st0_tag = FPU_Special(st0_ptr);
1337         if (st1_tag == TAG_Special)
1338                 st1_tag = FPU_Special(st1_ptr);
1339
1340         if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
1341             || ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
1342             || ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
1343                 if (denormal_operand() < 0)
1344                         return;
1345
1346                 goto valid_yl2xp1;
1347         } else if ((st0_tag == TAG_Empty) | (st1_tag == TAG_Empty)) {
1348                 FPU_stack_underflow_pop(1);
1349                 return;
1350         } else if (st0_tag == TAG_Zero) {
1351                 switch (st1_tag) {
1352                 case TW_Denormal:
1353                         if (denormal_operand() < 0)
1354                                 return;
1355                         /* fall through */
1356                 case TAG_Zero:
1357                 case TAG_Valid:
1358                         setsign(st0_ptr, getsign(st0_ptr) ^ getsign(st1_ptr));
1359                         FPU_copy_to_reg1(st0_ptr, st0_tag);
1360                         break;
1361
1362                 case TW_Infinity:
1363                         /* Infinity*log(1) */
1364                         if (arith_invalid(1) < 0)
1365                                 return;
1366                         break;
1367
1368                 case TW_NaN:
1369                         if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1370                                 return;
1371                         break;
1372
1373                 default:
1374 #ifdef PARANOID
1375                         EXCEPTION(EX_INTERNAL | 0x116);
1376                         return;
1377 #endif /* PARANOID */
1378                         break;
1379                 }
1380         } else if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
1381                 switch (st1_tag) {
1382                 case TAG_Zero:
1383                         if (signnegative(st0_ptr)) {
1384                                 if (exponent(st0_ptr) >= 0) {
1385                                         /* st(0) holds <= -1.0 */
1386 #ifdef PECULIAR_486             /* Stupid 80486 doesn't worry about log(negative). */
1387                                         changesign(st1_ptr);
1388 #else
1389                                         if (arith_invalid(1) < 0)
1390                                                 return;
1391 #endif /* PECULIAR_486 */
1392                                 } else if ((st0_tag == TW_Denormal)
1393                                            && (denormal_operand() < 0))
1394                                         return;
1395                                 else
1396                                         changesign(st1_ptr);
1397                         } else if ((st0_tag == TW_Denormal)
1398                                    && (denormal_operand() < 0))
1399                                 return;
1400                         break;
1401
1402                 case TW_Infinity:
1403                         if (signnegative(st0_ptr)) {
1404                                 if ((exponent(st0_ptr) >= 0) &&
1405                                     !((st0_ptr->sigh == 0x80000000) &&
1406                                       (st0_ptr->sigl == 0))) {
1407                                         /* st(0) holds < -1.0 */
1408 #ifdef PECULIAR_486             /* Stupid 80486 doesn't worry about log(negative). */
1409                                         changesign(st1_ptr);
1410 #else
1411                                         if (arith_invalid(1) < 0)
1412                                                 return;
1413 #endif /* PECULIAR_486 */
1414                                 } else if ((st0_tag == TW_Denormal)
1415                                            && (denormal_operand() < 0))
1416                                         return;
1417                                 else
1418                                         changesign(st1_ptr);
1419                         } else if ((st0_tag == TW_Denormal)
1420                                    && (denormal_operand() < 0))
1421                                 return;
1422                         break;
1423
1424                 case TW_NaN:
1425                         if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1426                                 return;
1427                 }
1428
1429         } else if (st0_tag == TW_NaN) {
1430                 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1431                         return;
1432         } else if (st0_tag == TW_Infinity) {
1433                 if (st1_tag == TW_NaN) {
1434                         if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1435                                 return;
1436                 } else if (signnegative(st0_ptr)) {
1437 #ifndef PECULIAR_486
1438                         /* This should have higher priority than denormals, but... */
1439                         if (arith_invalid(1) < 0)       /* log(-infinity) */
1440                                 return;
1441 #endif /* PECULIAR_486 */
1442                         if ((st1_tag == TW_Denormal)
1443                             && (denormal_operand() < 0))
1444                                 return;
1445 #ifdef PECULIAR_486
1446                         /* Denormal operands actually get higher priority */
1447                         if (arith_invalid(1) < 0)       /* log(-infinity) */
1448                                 return;
1449 #endif /* PECULIAR_486 */
1450                 } else if (st1_tag == TAG_Zero) {
1451                         /* log(infinity) */
1452                         if (arith_invalid(1) < 0)
1453                                 return;
1454                 }
1455
1456                 /* st(1) must be valid here. */
1457
1458                 else if ((st1_tag == TW_Denormal) && (denormal_operand() < 0))
1459                         return;
1460
1461                 /* The Manual says that log(Infinity) is invalid, but a real
1462                    80486 sensibly says that it is o.k. */
1463                 else {
1464                         u_char sign = getsign(st1_ptr);
1465                         FPU_copy_to_reg1(&CONST_INF, TAG_Special);
1466                         setsign(st1_ptr, sign);
1467                 }
1468         }
1469 #ifdef PARANOID
1470         else {
1471                 EXCEPTION(EX_INTERNAL | 0x117);
1472                 return;
1473         }
1474 #endif /* PARANOID */
1475
1476         FPU_pop();
1477         return;
1478
1479 }
1480
1481 static void fscale(FPU_REG *st0_ptr, u_char st0_tag)
1482 {
1483         FPU_REG *st1_ptr = &st(1);
1484         u_char st1_tag = FPU_gettagi(1);
1485         int old_cw = control_word;
1486         u_char sign = getsign(st0_ptr);
1487
1488         clear_C1();
1489         if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1490                 long scale;
1491                 FPU_REG tmp;
1492
1493                 /* Convert register for internal use. */
1494                 setexponent16(st0_ptr, exponent(st0_ptr));
1495
1496               valid_scale:
1497
1498                 if (exponent(st1_ptr) > 30) {
1499                         /* 2^31 is far too large, would require 2^(2^30) or 2^(-2^30) */
1500
1501                         if (signpositive(st1_ptr)) {
1502                                 EXCEPTION(EX_Overflow);
1503                                 FPU_copy_to_reg0(&CONST_INF, TAG_Special);
1504                         } else {
1505                                 EXCEPTION(EX_Underflow);
1506                                 FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
1507                         }
1508                         setsign(st0_ptr, sign);
1509                         return;
1510                 }
1511
1512                 control_word &= ~CW_RC;
1513                 control_word |= RC_CHOP;
1514                 reg_copy(st1_ptr, &tmp);
1515                 FPU_round_to_int(&tmp, st1_tag);        /* This can never overflow here */
1516                 control_word = old_cw;
1517                 scale = signnegative(st1_ptr) ? -tmp.sigl : tmp.sigl;
1518                 scale += exponent16(st0_ptr);
1519
1520                 setexponent16(st0_ptr, scale);
1521
1522                 /* Use FPU_round() to properly detect under/overflow etc */
1523                 FPU_round(st0_ptr, 0, 0, control_word, sign);
1524
1525                 return;
1526         }
1527
1528         if (st0_tag == TAG_Special)
1529                 st0_tag = FPU_Special(st0_ptr);
1530         if (st1_tag == TAG_Special)
1531                 st1_tag = FPU_Special(st1_ptr);
1532
1533         if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
1534                 switch (st1_tag) {
1535                 case TAG_Valid:
1536                         /* st(0) must be a denormal */
1537                         if ((st0_tag == TW_Denormal)
1538                             && (denormal_operand() < 0))
1539                                 return;
1540
1541                         FPU_to_exp16(st0_ptr, st0_ptr); /* Will not be left on stack */
1542                         goto valid_scale;
1543
1544                 case TAG_Zero:
1545                         if (st0_tag == TW_Denormal)
1546                                 denormal_operand();
1547                         return;
1548
1549                 case TW_Denormal:
1550                         denormal_operand();
1551                         return;
1552
1553                 case TW_Infinity:
1554                         if ((st0_tag == TW_Denormal)
1555                             && (denormal_operand() < 0))
1556                                 return;
1557
1558                         if (signpositive(st1_ptr))
1559                                 FPU_copy_to_reg0(&CONST_INF, TAG_Special);
1560                         else
1561                                 FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
1562                         setsign(st0_ptr, sign);
1563                         return;
1564
1565                 case TW_NaN:
1566                         real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1567                         return;
1568                 }
1569         } else if (st0_tag == TAG_Zero) {
1570                 switch (st1_tag) {
1571                 case TAG_Valid:
1572                 case TAG_Zero:
1573                         return;
1574
1575                 case TW_Denormal:
1576                         denormal_operand();
1577                         return;
1578
1579                 case TW_Infinity:
1580                         if (signpositive(st1_ptr))
1581                                 arith_invalid(0);       /* Zero scaled by +Infinity */
1582                         return;
1583
1584                 case TW_NaN:
1585                         real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1586                         return;
1587                 }
1588         } else if (st0_tag == TW_Infinity) {
1589                 switch (st1_tag) {
1590                 case TAG_Valid:
1591                 case TAG_Zero:
1592                         return;
1593
1594                 case TW_Denormal:
1595                         denormal_operand();
1596                         return;
1597
1598                 case TW_Infinity:
1599                         if (signnegative(st1_ptr))
1600                                 arith_invalid(0);       /* Infinity scaled by -Infinity */
1601                         return;
1602
1603                 case TW_NaN:
1604                         real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1605                         return;
1606                 }
1607         } else if (st0_tag == TW_NaN) {
1608                 if (st1_tag != TAG_Empty) {
1609                         real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1610                         return;
1611                 }
1612         }
1613 #ifdef PARANOID
1614         if (!((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty))) {
1615                 EXCEPTION(EX_INTERNAL | 0x115);
1616                 return;
1617         }
1618 #endif
1619
1620         /* At least one of st(0), st(1) must be empty */
1621         FPU_stack_underflow();
1622
1623 }
1624
1625 /*---------------------------------------------------------------------------*/
1626
1627 static FUNC_ST0 const trig_table_a[] = {
1628         f2xm1, fyl2x, fptan, fpatan,
1629         fxtract, fprem1, (FUNC_ST0) fdecstp, (FUNC_ST0) fincstp
1630 };
1631
1632 void FPU_triga(void)
1633 {
1634         (trig_table_a[FPU_rm]) (&st(0), FPU_gettag0());
1635 }
1636
1637 static FUNC_ST0 const trig_table_b[] = {
1638         fprem, fyl2xp1, fsqrt_, fsincos, frndint_, fscale, (FUNC_ST0) fsin, fcos
1639 };
1640
1641 void FPU_trigb(void)
1642 {
1643         (trig_table_b[FPU_rm]) (&st(0), FPU_gettag0());
1644 }