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