开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > Linux/Unix编程 > 查看源码
fpu_trig.c
资源名称：linux-2.4.20.tar.gz　[点击查看]
上传用户：jlfgdled
上传日期：2013-04-10
资源大小：33168k
文件大小：41k
源码类别：

Linux/Unix编程

开发平台：

Unix_Linux

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