op-2.h
上传用户:jlfgdled
上传日期:2013-04-10
资源大小:33168k
文件大小:13k
源码类别:

Linux/Unix编程

开发平台:

Unix_Linux

  1. /*
  2.  * BK Id: SCCS/s.op-2.h 1.5 05/17/01 18:14:23 cort
  3.  */
  4. /*
  5.  * Basic two-word fraction declaration and manipulation.
  6.  */
  7. #define _FP_FRAC_DECL_2(X) _FP_W_TYPE X##_f0, X##_f1
  8. #define _FP_FRAC_COPY_2(D,S) (D##_f0 = S##_f0, D##_f1 = S##_f1)
  9. #define _FP_FRAC_SET_2(X,I) __FP_FRAC_SET_2(X, I)
  10. #define _FP_FRAC_HIGH_2(X) (X##_f1)
  11. #define _FP_FRAC_LOW_2(X) (X##_f0)
  12. #define _FP_FRAC_WORD_2(X,w) (X##_f##w)
  13. #define _FP_FRAC_SLL_2(X,N)
  14.   do {
  15.     if ((N) < _FP_W_TYPE_SIZE)
  16.       {
  17.         if (__builtin_constant_p(N) && (N) == 1) 
  18.           {
  19.             X##_f1 = X##_f1 + X##_f1 + (((_FP_WS_TYPE)(X##_f0)) < 0);
  20.             X##_f0 += X##_f0;
  21.           }
  22.         else
  23.           {
  24.     X##_f1 = X##_f1 << (N) | X##_f0 >> (_FP_W_TYPE_SIZE - (N));
  25.     X##_f0 <<= (N);
  26.   }
  27.       }
  28.     else
  29.       {
  30. X##_f1 = X##_f0 << ((N) - _FP_W_TYPE_SIZE);
  31. X##_f0 = 0;
  32.       }
  33.   } while (0)
  34. #define _FP_FRAC_SRL_2(X,N)
  35.   do {
  36.     if ((N) < _FP_W_TYPE_SIZE)
  37.       {
  38. X##_f0 = X##_f0 >> (N) | X##_f1 << (_FP_W_TYPE_SIZE - (N));
  39. X##_f1 >>= (N);
  40.       }
  41.     else
  42.       {
  43. X##_f0 = X##_f1 >> ((N) - _FP_W_TYPE_SIZE);
  44. X##_f1 = 0;
  45.       }
  46.   } while (0)
  47. /* Right shift with sticky-lsb.  */
  48. #define _FP_FRAC_SRS_2(X,N,sz)
  49.   do {
  50.     if ((N) < _FP_W_TYPE_SIZE)
  51.       {
  52. X##_f0 = (X##_f1 << (_FP_W_TYPE_SIZE - (N)) | X##_f0 >> (N) |
  53.   (__builtin_constant_p(N) && (N) == 1
  54.    ? X##_f0 & 1
  55.    : (X##_f0 << (_FP_W_TYPE_SIZE - (N))) != 0));
  56. X##_f1 >>= (N);
  57.       }
  58.     else
  59.       {
  60. X##_f0 = (X##_f1 >> ((N) - _FP_W_TYPE_SIZE) |
  61.           (((X##_f1 << (sz - (N))) | X##_f0) != 0));
  62. X##_f1 = 0;
  63.       }
  64.   } while (0)
  65. #define _FP_FRAC_ADDI_2(X,I) 
  66.   __FP_FRAC_ADDI_2(X##_f1, X##_f0, I)
  67. #define _FP_FRAC_ADD_2(R,X,Y) 
  68.   __FP_FRAC_ADD_2(R##_f1, R##_f0, X##_f1, X##_f0, Y##_f1, Y##_f0)
  69. #define _FP_FRAC_SUB_2(R,X,Y) 
  70.   __FP_FRAC_SUB_2(R##_f1, R##_f0, X##_f1, X##_f0, Y##_f1, Y##_f0)
  71. #define _FP_FRAC_CLZ_2(R,X)
  72.   do {
  73.     if (X##_f1)
  74.       __FP_CLZ(R,X##_f1);
  75.     else 
  76.     {
  77.       __FP_CLZ(R,X##_f0);
  78.       R += _FP_W_TYPE_SIZE;
  79.     }
  80.   } while(0)
  81. /* Predicates */
  82. #define _FP_FRAC_NEGP_2(X) ((_FP_WS_TYPE)X##_f1 < 0)
  83. #define _FP_FRAC_ZEROP_2(X) ((X##_f1 | X##_f0) == 0)
  84. #define _FP_FRAC_OVERP_2(fs,X) (X##_f1 & _FP_OVERFLOW_##fs)
  85. #define _FP_FRAC_EQ_2(X, Y) (X##_f1 == Y##_f1 && X##_f0 == Y##_f0)
  86. #define _FP_FRAC_GT_2(X, Y)
  87.   ((X##_f1 > Y##_f1) || (X##_f1 == Y##_f1 && X##_f0 > Y##_f0))
  88. #define _FP_FRAC_GE_2(X, Y)
  89.   ((X##_f1 > Y##_f1) || (X##_f1 == Y##_f1 && X##_f0 >= Y##_f0))
  90. #define _FP_ZEROFRAC_2 0, 0
  91. #define _FP_MINFRAC_2 0, 1
  92. /*
  93.  * Internals 
  94.  */
  95. #define __FP_FRAC_SET_2(X,I1,I0) (X##_f0 = I0, X##_f1 = I1)
  96. #define __FP_CLZ_2(R, xh, xl)
  97.   do {
  98.     if (xh)
  99.       __FP_CLZ(R,xl);
  100.     else 
  101.     {
  102.       __FP_CLZ(R,xl);
  103.       R += _FP_W_TYPE_SIZE;
  104.     }
  105.   } while(0)
  106. #if 0
  107. #ifndef __FP_FRAC_ADDI_2
  108. #define __FP_FRAC_ADDI_2(xh, xl, i) 
  109.   (xh += ((xl += i) < i))
  110. #endif
  111. #ifndef __FP_FRAC_ADD_2
  112. #define __FP_FRAC_ADD_2(rh, rl, xh, xl, yh, yl) 
  113.   (rh = xh + yh + ((rl = xl + yl) < xl))
  114. #endif
  115. #ifndef __FP_FRAC_SUB_2
  116. #define __FP_FRAC_SUB_2(rh, rl, xh, xl, yh, yl) 
  117.   (rh = xh - yh - ((rl = xl - yl) > xl))
  118. #endif
  119. #else
  120. #undef __FP_FRAC_ADDI_2
  121. #define __FP_FRAC_ADDI_2(xh, xl, i) add_ssaaaa(xh, xl, xh, xl, 0, i)
  122. #undef __FP_FRAC_ADD_2
  123. #define __FP_FRAC_ADD_2 add_ssaaaa
  124. #undef __FP_FRAC_SUB_2
  125. #define __FP_FRAC_SUB_2 sub_ddmmss
  126. #endif
  127. /*
  128.  * Unpack the raw bits of a native fp value.  Do not classify or
  129.  * normalize the data.
  130.  */
  131. #define _FP_UNPACK_RAW_2(fs, X, val)
  132.   do {
  133.     union _FP_UNION_##fs _flo; _flo.flt = (val);
  134.     X##_f0 = _flo.bits.frac0;
  135.     X##_f1 = _flo.bits.frac1;
  136.     X##_e  = _flo.bits.exp;
  137.     X##_s  = _flo.bits.sign;
  138.   } while (0)
  139. /*
  140.  * Repack the raw bits of a native fp value.
  141.  */
  142. #define _FP_PACK_RAW_2(fs, val, X)
  143.   do {
  144.     union _FP_UNION_##fs _flo;
  145.     _flo.bits.frac0 = X##_f0;
  146.     _flo.bits.frac1 = X##_f1;
  147.     _flo.bits.exp   = X##_e;
  148.     _flo.bits.sign  = X##_s;
  149.     (val) = _flo.flt;
  150.   } while (0)
  151. /*
  152.  * Multiplication algorithms:
  153.  */
  154. /* Given a 1W * 1W => 2W primitive, do the extended multiplication.  */
  155. #define _FP_MUL_MEAT_2_wide(fs, R, X, Y, doit)
  156.   do {
  157.     _FP_FRAC_DECL_4(_z); _FP_FRAC_DECL_2(_b); _FP_FRAC_DECL_2(_c);
  158.     doit(_FP_FRAC_WORD_4(_z,1), _FP_FRAC_WORD_4(_z,0), X##_f0, Y##_f0); 
  159.     doit(_b_f1, _b_f0, X##_f0, Y##_f1);
  160.     doit(_c_f1, _c_f0, X##_f1, Y##_f0);
  161.     doit(_FP_FRAC_WORD_4(_z,3), _FP_FRAC_WORD_4(_z,2), X##_f1, Y##_f1); 
  162.     __FP_FRAC_ADD_4(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),
  163.     _FP_FRAC_WORD_4(_z,1),_FP_FRAC_WORD_4(_z,0),
  164.     0, _b_f1, _b_f0, 0,
  165.     _FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),
  166.     _FP_FRAC_WORD_4(_z,1),_FP_FRAC_WORD_4(_z,0));
  167.     __FP_FRAC_ADD_4(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),
  168.     _FP_FRAC_WORD_4(_z,1),_FP_FRAC_WORD_4(_z,0),
  169.     0, _c_f1, _c_f0, 0,
  170.     _FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),
  171.     _FP_FRAC_WORD_4(_z,1),_FP_FRAC_WORD_4(_z,0));
  172.     /* Normalize since we know where the msb of the multiplicands
  173.        were (bit B), we know that the msb of the of the product is
  174.        at either 2B or 2B-1.  */
  175.     _FP_FRAC_SRS_4(_z, _FP_WFRACBITS_##fs-1, 2*_FP_WFRACBITS_##fs);
  176.     R##_f0 = _FP_FRAC_WORD_4(_z,0);
  177.     R##_f1 = _FP_FRAC_WORD_4(_z,1);
  178.   } while (0)
  179. /* This next macro appears to be totally broken. Fortunately nowhere
  180.  * seems to use it :-> The problem is that we define _z[4] but
  181.  * then use it in _FP_FRAC_SRS_4, which will attempt to access
  182.  * _z_f[n] which will cause an error. The fix probably involves 
  183.  * declaring it with _FP_FRAC_DECL_4, see previous macro. -- PMM 02/1998 
  184.  */
  185. #define _FP_MUL_MEAT_2_gmp(fs, R, X, Y)
  186.   do {
  187.     _FP_W_TYPE _x[2], _y[2], _z[4];
  188.     _x[0] = X##_f0; _x[1] = X##_f1;
  189.     _y[0] = Y##_f0; _y[1] = Y##_f1;
  190.     mpn_mul_n(_z, _x, _y, 2);
  191.     /* Normalize since we know where the msb of the multiplicands
  192.        were (bit B), we know that the msb of the of the product is
  193.        at either 2B or 2B-1.  */
  194.     _FP_FRAC_SRS_4(_z, _FP_WFRACBITS##_fs-1, 2*_FP_WFRACBITS_##fs);
  195.     R##_f0 = _z[0];
  196.     R##_f1 = _z[1];
  197.   } while (0)
  198. /*
  199.  * Division algorithms:
  200.  * This seems to be giving me difficulties -- PMM 
  201.  * Look, NetBSD seems to be able to comment algorithms. Can't you?
  202.  * I've thrown printks at the problem.
  203.  * This now appears to work, but I still don't really know why.
  204.  * Also, I don't think the result is properly normalised...
  205.  */
  206. #define _FP_DIV_MEAT_2_udiv_64(fs, R, X, Y)
  207.   do {
  208.     extern void _fp_udivmodti4(_FP_W_TYPE q[2], _FP_W_TYPE r[2],
  209.        _FP_W_TYPE n1, _FP_W_TYPE n0,
  210.        _FP_W_TYPE d1, _FP_W_TYPE d0);
  211.     _FP_W_TYPE _n_f3, _n_f2, _n_f1, _n_f0, _r_f1, _r_f0;
  212.     _FP_W_TYPE _q_f1, _q_f0, _m_f1, _m_f0;
  213.     _FP_W_TYPE _rmem[2], _qmem[2];
  214.     /* I think this check is to ensure that the result is normalised.   
  215.      * Assuming X,Y normalised (ie in [1.0,2.0)) X/Y will be in         
  216.      * [0.5,2.0). Furthermore, it will be less than 1.0 iff X < Y.      
  217.      * In this case we tweak things. (this is based on comments in      
  218.      * the NetBSD FPU emulation code. )                                 
  219.      * We know X,Y are normalised because we ensure this as part of     
  220.      * the unpacking process. -- PMM                                    
  221.      */
  222.     if (_FP_FRAC_GT_2(X, Y))
  223.       {
  224. /* R##_e++; */
  225. _n_f3 = X##_f1 >> 1;
  226. _n_f2 = X##_f1 << (_FP_W_TYPE_SIZE - 1) | X##_f0 >> 1;
  227. _n_f1 = X##_f0 << (_FP_W_TYPE_SIZE - 1);
  228. _n_f0 = 0;
  229.       }
  230.     else
  231.       {
  232. R##_e--;
  233. _n_f3 = X##_f1;
  234. _n_f2 = X##_f0;
  235. _n_f1 = _n_f0 = 0;
  236.       }
  237.     /* Normalize, i.e. make the most significant bit of the 
  238.        denominator set.  CHANGED: - 1 to nothing -- PMM */
  239.     _FP_FRAC_SLL_2(Y, _FP_WFRACXBITS_##fs /* -1 */);
  240.     /* Do the 256/128 bit division given the 128-bit _fp_udivmodtf4 
  241.        primitive snagged from libgcc2.c.  */
  242.     _fp_udivmodti4(_qmem, _rmem, _n_f3, _n_f2, 0, Y##_f1);
  243.     _q_f1 = _qmem[0];
  244.     umul_ppmm(_m_f1, _m_f0, _q_f1, Y##_f0);
  245.     _r_f1 = _rmem[0];
  246.     _r_f0 = _n_f1;
  247.     if (_FP_FRAC_GT_2(_m, _r))
  248.       {
  249. _q_f1--;
  250. _FP_FRAC_ADD_2(_r, _r, Y);
  251. if (_FP_FRAC_GE_2(_r, Y) && _FP_FRAC_GT_2(_m, _r))
  252.   {
  253.     _q_f1--;
  254.     _FP_FRAC_ADD_2(_r, _r, Y);
  255.   }
  256.       }
  257.     _FP_FRAC_SUB_2(_r, _r, _m);
  258.     _fp_udivmodti4(_qmem, _rmem, _r_f1, _r_f0, 0, Y##_f1);
  259.     _q_f0 = _qmem[0];
  260.     umul_ppmm(_m_f1, _m_f0, _q_f0, Y##_f0);
  261.     _r_f1 = _rmem[0];
  262.     _r_f0 = _n_f0;
  263.     if (_FP_FRAC_GT_2(_m, _r))
  264.       {
  265. _q_f0--;
  266. _FP_FRAC_ADD_2(_r, _r, Y);
  267. if (_FP_FRAC_GE_2(_r, Y) && _FP_FRAC_GT_2(_m, _r))
  268.   {
  269.     _q_f0--;
  270.     _FP_FRAC_ADD_2(_r, _r, Y);
  271.   }
  272.       }
  273.     _FP_FRAC_SUB_2(_r, _r, _m);
  274.     R##_f1 = _q_f1;
  275.     R##_f0 = _q_f0 | ((_r_f1 | _r_f0) != 0);
  276.     /* adjust so answer is normalized again. I'm not sure what the 
  277.      * final sz param should be. In practice it's never used since      
  278.      * N is 1 which is always going to be < _FP_W_TYPE_SIZE...
  279.      */
  280.     /* _FP_FRAC_SRS_2(R,1,_FP_WFRACBITS_##fs); */
  281.   } while (0)
  282. #define _FP_DIV_MEAT_2_gmp(fs, R, X, Y)
  283.   do {
  284.     _FP_W_TYPE _x[4], _y[2], _z[4];
  285.     _y[0] = Y##_f0; _y[1] = Y##_f1;
  286.     _x[0] = _x[3] = 0;
  287.     if (_FP_FRAC_GT_2(X, Y))
  288.       {
  289. R##_e++;
  290. _x[1] = (X##_f0 << (_FP_WFRACBITS-1 - _FP_W_TYPE_SIZE) |
  291.  X##_f1 >> (_FP_W_TYPE_SIZE -
  292.     (_FP_WFRACBITS-1 - _FP_W_TYPE_SIZE)));
  293. _x[2] = X##_f1 << (_FP_WFRACBITS-1 - _FP_W_TYPE_SIZE);
  294.       }
  295.     else
  296.       {
  297. _x[1] = (X##_f0 << (_FP_WFRACBITS - _FP_W_TYPE_SIZE) |
  298.  X##_f1 >> (_FP_W_TYPE_SIZE -
  299.     (_FP_WFRACBITS - _FP_W_TYPE_SIZE)));
  300. _x[2] = X##_f1 << (_FP_WFRACBITS - _FP_W_TYPE_SIZE);
  301.       }
  302.     (void) mpn_divrem (_z, 0, _x, 4, _y, 2);
  303.     R##_f1 = _z[1];
  304.     R##_f0 = _z[0] | ((_x[0] | _x[1]) != 0);
  305.   } while (0)
  306. /*
  307.  * Square root algorithms:
  308.  * We have just one right now, maybe Newton approximation
  309.  * should be added for those machines where division is fast.
  310.  */
  311.  
  312. #define _FP_SQRT_MEAT_2(R, S, T, X, q)
  313.   do {
  314.     while (q)
  315.       {
  316.         T##_f1 = S##_f1 + q;
  317.         if (T##_f1 <= X##_f1)
  318.           {
  319.             S##_f1 = T##_f1 + q;
  320.             X##_f1 -= T##_f1;
  321.             R##_f1 += q;
  322.           }
  323.         _FP_FRAC_SLL_2(X, 1);
  324.         q >>= 1;
  325.       }
  326.     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);
  327.     while (q)
  328.       {
  329.         T##_f0 = S##_f0 + q;
  330.         T##_f1 = S##_f1;
  331.         if (T##_f1 < X##_f1 || 
  332.             (T##_f1 == X##_f1 && T##_f0 < X##_f0))
  333.           {
  334.             S##_f0 = T##_f0 + q;
  335.             if (((_FP_WS_TYPE)T##_f0) < 0 &&
  336.                 ((_FP_WS_TYPE)S##_f0) >= 0)
  337.               S##_f1++;
  338.             _FP_FRAC_SUB_2(X, X, T);
  339.             R##_f0 += q;
  340.           }
  341.         _FP_FRAC_SLL_2(X, 1);
  342.         q >>= 1;
  343.       }
  344.   } while (0)
  345. /*
  346.  * Assembly/disassembly for converting to/from integral types.  
  347.  * No shifting or overflow handled here.
  348.  */
  349. #define _FP_FRAC_ASSEMBLE_2(r, X, rsize)
  350.   do {
  351.     if (rsize <= _FP_W_TYPE_SIZE)
  352.       r = X##_f0;
  353.     else
  354.       {
  355. r = X##_f1;
  356. r <<= _FP_W_TYPE_SIZE;
  357. r += X##_f0;
  358.       }
  359.   } while (0)
  360. #define _FP_FRAC_DISASSEMBLE_2(X, r, rsize)
  361.   do {
  362.     X##_f0 = r;
  363.     X##_f1 = (rsize <= _FP_W_TYPE_SIZE ? 0 : r >> _FP_W_TYPE_SIZE);
  364.   } while (0)
  365. /*
  366.  * Convert FP values between word sizes
  367.  */
  368. #define _FP_FRAC_CONV_1_2(dfs, sfs, D, S)
  369.   do {
  370.     _FP_FRAC_SRS_2(S, (_FP_WFRACBITS_##sfs - _FP_WFRACBITS_##dfs),
  371.    _FP_WFRACBITS_##sfs);
  372.     D##_f = S##_f0;
  373.   } while (0)
  374. #define _FP_FRAC_CONV_2_1(dfs, sfs, D, S)
  375.   do {
  376.     D##_f0 = S##_f;
  377.     D##_f1 = 0;
  378.     _FP_FRAC_SLL_2(D, (_FP_WFRACBITS_##dfs - _FP_WFRACBITS_##sfs));
  379.   } while (0)