home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk1.iso / altsrc / articles / 11179 < prev    next >
Text File  |  1994-08-26  |  5KB  |  164 lines
  1. Newsgroups: comp.sys.powerpc,comp.arch.arithmetic,gnu.misc.discuss,comp.programming,sci.math,alt.sources
  2. Path: wupost!gumby!newsxfer.itd.umich.edu!isclient.merit.edu!msuinfo!agate!howland.reston.ans.net!cs.utexas.edu!rutgers!att-out!nntpa!babel.ho.att.com!joe
  3. From: joe@sanskrit.ho.att.com (Joe Orost)
  4. Subject: ANSWER: algorithm to perform 64-bit / 32-bit signed & unsigned divide
  5. Message-ID: <joe.777905986@babel.ho.att.com>
  6. Keywords: division
  7. Sender: news@nntpa.cb.att.com (Netnews Administration)
  8. Nntp-Posting-Host: sanskrit.ho.att.com
  9. Organization: AT&T
  10. References: <joe.777262055@babel.ho.att.com>
  11. Date: Fri, 26 Aug 1994 12:59:46 GMT
  12. Lines: 149
  13. Xref: wupost comp.sys.powerpc:27845 comp.arch.arithmetic:571 gnu.misc.discuss:18291 comp.programming:12089 sci.math:75673 alt.sources:11179
  14.  
  15. /*
  16.  * Unsigned divide 64-bit / 32-bit
  17.  *
  18.  * Input:
  19.  *     uint D    32-bit divisor
  20.  *     uint H    upper 32-bits of dividend
  21.  *     uint L    lower 32-bits of dividend
  22.  *     uint *RP   pointer to 32 bit remainder
  23.  *     uint *QP   pointer to 32 bit quotient
  24.  *
  25.  * Output:
  26.  *     int       overflow flag
  27.  *     *RP       32-bit remainder
  28.  *     *QP       32-bit quotient
  29.  *
  30.  * Notes:
  31.  *     1) If quotient does not fit in 32-bits then 1 is returned
  32.  *
  33.  *     2) The D2INT macro truncates a double to an integer.
  34.  *        It is coded for BIG ENDIAN.
  35.  *        The rounding mode must be set to round towards zero.
  36.  *
  37.  * Authors: Mauricio Breternitz Jr.     (Motorola)
  38.  *          Arturo Martin-de-Nicolas    (IBM)
  39.  *          Joseph M. Orost             (AT&T Bell Labs)
  40.  *
  41.  * Algorithm:
  42.  *      1) Verify that the quotient will fit in 32-bits.
  43.  *      2) Compute the reciprocal of the divisor using
  44.  *         the floating point unit.
  45.  *      3) Compute an approximation of the quotient by
  46.  *         multiplying each half of the dividend by
  47.  *         the reciprocal.
  48.  *      4) Truncate the approximate quotient to an
  49.  *         integer.  Due to rounding, this value may be
  50.  *         smaller than the actual quotient by one.
  51.  *      5) Compute the remainder based on the approximate
  52.  *         quotient.
  53.  *      6) If necessary, increment quotient and adjust
  54.  *         remainder.  There are two cases of an incorrect
  55.  *         quotient:
  56.  *              a) The computed remainder is greater than
  57.  *                 or equal to the divisor.
  58.  *              b) The computed remainder is smaller than
  59.  *                 an approximate remainder computed by
  60.  *                 multiplying the fractional part of the
  61.  *                 quotient with the divisor.
  62.  */
  63. /************************************************************************/
  64. /*  WE PROVIDE THIS SOURCE CODE EXAMPLE "AS IS," WITHOUT WARRANTY OF    */
  65. /*  ANY KIND EITHER EXPRESSED OR IMPLIED INCLUDING BUT NOT LIMITED TO   */
  66. /*  THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A         */
  67. /*  PARTICULAR PURPOSE.  The entire risk as to the quality and          */
  68. /*  performance of this example is with you.                            */
  69. /************************************************************************/
  70.  
  71. /*
  72.  *
  73.  *      Values for the IEEE Rounding Modes (ANSI Encoding)
  74.  *
  75.  *      RZ = Round toward zero
  76.  *      RN = Round toward nearest (default)
  77.  *      RP = Round toward plus infinity
  78.  *      RM = Round toward minus infinity
  79.  *
  80.  */
  81. #define FP_RND_RZ       0
  82. #define FP_RND_RN       1
  83. #define FP_RND_RP       2
  84. #define FP_RND_RM       3
  85.  
  86.  
  87. #define TWO32 ((double)(1<<16)*(double)(1<<16))
  88. #define TWO52 ((double)(1<<26)*(double)(1<<26))
  89.  
  90. #define uint unsigned int
  91.  
  92. #define D2INT(DOUBLE,INT) { double convert = (DOUBLE) + TWO52; \
  93.                             INT = ((uint *)(&convert))[1]; }
  94.  
  95.  
  96. int divu_64_32( uint D, uint H, uint L, uint *RP, uint *QP )
  97. {
  98.     uint Q;
  99.  
  100.     /* Determine if quotient fits in 32-bits */
  101.     if (H >= D)
  102.         return 1;         /* Overflow */
  103.     else {
  104.         unsigned short save = fp_swap_rnd( FP_RND_RZ );
  105.         double rd = 1 / (double)D;
  106.         double Q_f = rd * (double)H * TWO32 + rd * (double)L;
  107.         D2INT( Q_f, Q );
  108.         {
  109.             uint R_aprox;
  110.             uint R_qd = L - Q * D;
  111.             if (R_qd >= D)
  112.                 Q++, *RP = R_qd-D;
  113.             else {
  114.                 D2INT((Q_f - (double)Q) * (double)D, R_aprox);
  115.                 if (R_aprox > R_qd)
  116.                     Q++, *RP = R_qd-D;
  117.                 else
  118.                     *RP = R_qd;
  119.             }
  120.         }
  121.         (void) fp_swap_rnd( save );
  122.     }
  123.  
  124.     *QP = Q;
  125.     return 0;
  126. }
  127.  
  128. /* The rest of this is written by Joe Orost (joe@babel.ho.att.com) */
  129.  
  130. #define sign(a) ((a) < 0)
  131. #define abs(a)  (sign(a) ? -(a) : (a))
  132.  
  133. int
  134. divs_64_32(d, n_0, n_1, RP, QP)
  135. long d, n_0;
  136. unsigned long n_1;
  137. long *RP, *QP;
  138. {
  139.     register long s_n, s_d, ov;
  140.     unsigned long a_n_0, a_d;
  141.     register long long nn, a_nn;            /* Compile with gcc */
  142.  
  143.     nn = (((long long)(n_0)) << 32) + n_1;
  144.     s_n = sign(nn);
  145.     a_nn = abs(nn);
  146.     a_n_0 = (long)(a_nn >> 32);
  147.     n_1 = (long)(a_nn & 0xffffffff);
  148.     s_d = sign(d);
  149.     a_d = abs(d);
  150.     ov = divu_64_32(a_d, a_n_0, n_1, RP, QP);
  151.     if(ov)
  152.         return 1;
  153.     if(s_n != s_d) {
  154.         *QP = -(*QP);
  155.         if(*QP > 0)
  156.             return 1;
  157.     } else if(*QP < 0) {
  158.         return 1;
  159.     }
  160.     if(s_n)
  161.         *RP = -(*RP);
  162.     return 0;
  163. }
  164.