home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / gmp202.zip / mpn / generic / gcd.c < prev    next >
C/C++ Source or Header  |  1996-06-04  |  12KB  |  403 lines
  1. /* mpn/gcd.c: mpn_gcd for gcd of two odd integers.
  2.  
  3. Copyright (C) 1991, 1993, 1994, 1995, 1996 Free Software Foundation, Inc.
  4.  
  5. This file is part of the GNU MP Library.
  6.  
  7. The GNU MP Library is free software; you can redistribute it and/or modify
  8. it under the terms of the GNU Library General Public License as published by
  9. the Free Software Foundation; either version 2 of the License, or (at your
  10. option) any later version.
  11.  
  12. The GNU MP Library is distributed in the hope that it will be useful, but
  13. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
  15. License for more details.
  16.  
  17. You should have received a copy of the GNU Library General Public License
  18. along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
  19. the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  20. MA 02111-1307, USA. */
  21.  
  22. /* Integer greatest common divisor of two unsigned integers, using
  23.    the accelerated algorithm (see reference below).
  24.  
  25.    mp_size_t mpn_gcd (vp, vsize, up, usize).
  26.  
  27.    Preconditions [U = (up, usize) and V = (vp, vsize)]:
  28.  
  29.    1.  V is odd.
  30.    2.  numbits(U) >= numbits(V).
  31.  
  32.    Both U and V are destroyed by the operation.  The result is left at vp,
  33.    and its size is returned.
  34.  
  35.    Ken Weber (kweber@mat.ufrgs.br, kweber@mcs.kent.edu)
  36.  
  37.    Funding for this work has been partially provided by Conselho Nacional
  38.    de Desenvolvimento Cienti'fico e Tecnolo'gico (CNPq) do Brazil, Grant
  39.    301314194-2, and was done while I was a visiting reseacher in the Instituto
  40.    de Matema'tica at Universidade Federal do Rio Grande do Sul (UFRGS).
  41.  
  42.    Refer to
  43.     K. Weber, The accelerated integer GCD algorithm, ACM Transactions on
  44.     Mathematical Software, v. 21 (March), 1995, pp. 111-122.  */
  45.  
  46. #include "gmp.h"
  47. #include "gmp-impl.h"
  48. #include "longlong.h"
  49.  
  50. /* If MIN (usize, vsize) > ACCEL_THRESHOLD, then the accelerated algorithm is
  51.    used, otherwise the binary algorithm is used.  This may be adjusted for
  52.    different architectures.  */
  53. #ifndef ACCEL_THRESHOLD
  54. #define ACCEL_THRESHOLD 4
  55. #endif
  56.  
  57. /* When U and V differ in size by more than BMOD_THRESHOLD, the accelerated
  58.    algorithm reduces using the bmod operation.  Otherwise, the k-ary reduction
  59.    is used.  0 <= BMOD_THRESHOLD < BITS_PER_MP_LIMB.  */
  60. enum
  61.   {
  62.     BMOD_THRESHOLD = BITS_PER_MP_LIMB/2
  63.   };
  64.  
  65. #define SIGN_BIT  (~(~(mp_limb_t)0 >> 1))
  66.  
  67.  
  68. #define SWAP_LIMB(UL, VL) do{mp_limb_t __l=(UL);(UL)=(VL);(VL)=__l;}while(0)
  69. #define SWAP_PTR(UP, VP) do{mp_ptr __p=(UP);(UP)=(VP);(VP)=__p;}while(0)
  70. #define SWAP_SZ(US, VS) do{mp_size_t __s=(US);(US)=(VS);(VS)=__s;}while(0)
  71. #define SWAP_MPN(UP, US, VP, VS) do{SWAP_PTR(UP,VP);SWAP_SZ(US,VS);}while(0)
  72.  
  73. /* Use binary algorithm to compute V <-- GCD (V, U) for usize, vsize == 2.
  74.    Both U and V must be odd.  */
  75. static __gmp_inline mp_size_t
  76. #if __STDC__
  77. gcd_2 (mp_ptr vp, mp_srcptr up)
  78. #else
  79. gcd_2 (vp, up)
  80.      mp_ptr vp;
  81.      mp_srcptr up;
  82. #endif
  83. {
  84.   mp_limb_t u0, u1, v0, v1;
  85.   mp_size_t vsize;
  86.  
  87.   u0 = up[0], u1 = up[1], v0 = vp[0], v1 = vp[1];
  88.  
  89.   while (u1 != v1 && u0 != v0)
  90.     {
  91.       unsigned long int r;
  92.       if (u1 > v1)
  93.     {
  94.       u1 -= v1 + (u0 < v0), u0 -= v0;
  95.       count_trailing_zeros (r, u0);
  96.       u0 = u1 << (BITS_PER_MP_LIMB - r) | u0 >> r;
  97.       u1 >>= r;
  98.     }
  99.       else  /* u1 < v1.  */
  100.     {
  101.       v1 -= u1 + (v0 < u0), v0 -= u0;
  102.       count_trailing_zeros (r, v0);
  103.       v0 = v1 << (BITS_PER_MP_LIMB - r) | v0 >> r;
  104.       v1 >>= r;
  105.     }
  106.     }
  107.  
  108.   vp[0] = v0, vp[1] = v1, vsize = 1 + (v1 != 0);
  109.  
  110.   /* If U == V == GCD, done.  Otherwise, compute GCD (V, |U - V|).  */
  111.   if (u1 == v1 && u0 == v0)
  112.     return vsize;
  113.  
  114.   v0 = (u0 == v0) ? (u1 > v1) ? u1-v1 : v1-u1 : (u0 > v0) ? u0-v0 : v0-u0;
  115.   vp[0] = mpn_gcd_1 (vp, vsize, v0);
  116.  
  117.   return 1;
  118. }
  119.  
  120. /* The function find_a finds 0 < N < 2^BITS_PER_MP_LIMB such that there exists
  121.    0 < |D| < 2^BITS_PER_MP_LIMB, and N == D * C mod 2^(2*BITS_PER_MP_LIMB).
  122.    In the reference article, D was computed along with N, but it is better to
  123.    compute D separately as D <-- N / C mod 2^(BITS_PER_MP_LIMB + 1), treating
  124.    the result as a twos' complement signed integer.
  125.  
  126.    Initialize N1 to C mod 2^(2*BITS_PER_MP_LIMB).  According to the reference
  127.    article, N2 should be initialized to 2^(2*BITS_PER_MP_LIMB), but we use
  128.    2^(2*BITS_PER_MP_LIMB) - N1 to start the calculations within double
  129.    precision.  If N2 > N1 initially, the first iteration of the while loop
  130.    will swap them.  In all other situations, N1 >= N2 is maintained.  */
  131.  
  132. static __gmp_inline mp_limb_t
  133. #if __STDC__
  134. find_a (mp_srcptr cp)
  135. #else
  136. find_a (cp)
  137.      mp_srcptr cp;
  138. #endif
  139. {
  140.   unsigned long int leading_zero_bits = 0;
  141.  
  142.   mp_limb_t n1_l = cp[0];    /* N1 == n1_h * 2^BITS_PER_MP_LIMB + n1_l.  */
  143.   mp_limb_t n1_h = cp[1];
  144.  
  145.   mp_limb_t n2_l = -n1_l;    /* N2 == n2_h * 2^BITS_PER_MP_LIMB + n2_l.  */
  146.   mp_limb_t n2_h = ~n1_h;
  147.  
  148.   /* Main loop.  */
  149.   while (n2_h)            /* While N2 >= 2^BITS_PER_MP_LIMB.  */
  150.     {
  151.       /* N1 <-- N1 % N2.  */
  152.       if ((SIGN_BIT >> leading_zero_bits & n2_h) == 0)
  153.     {
  154.       unsigned long int i;
  155.       count_leading_zeros (i, n2_h);
  156.       i -= leading_zero_bits, leading_zero_bits += i;
  157.       n2_h = n2_h<<i | n2_l>>(BITS_PER_MP_LIMB - i), n2_l <<= i;
  158.       do
  159.         {
  160.           if (n1_h > n2_h || (n1_h == n2_h && n1_l >= n2_l))
  161.         n1_h -= n2_h + (n1_l < n2_l), n1_l -= n2_l;
  162.           n2_l = n2_l>>1 | n2_h<<(BITS_PER_MP_LIMB - 1), n2_h >>= 1;
  163.           i -= 1;
  164.         }
  165.       while (i);
  166.     }
  167.       if (n1_h > n2_h || (n1_h == n2_h && n1_l >= n2_l))
  168.     n1_h -= n2_h + (n1_l < n2_l), n1_l -= n2_l;
  169.  
  170.       SWAP_LIMB (n1_h, n2_h);
  171.       SWAP_LIMB (n1_l, n2_l);
  172.     }
  173.  
  174.   return n2_l;
  175. }
  176.  
  177. mp_size_t
  178. #if __STDC__
  179. mpn_gcd (mp_ptr gp, mp_ptr vp, mp_size_t vsize, mp_ptr up, mp_size_t usize)
  180. #else
  181. mpn_gcd (gp, vp, vsize, up, usize)
  182.      mp_ptr gp;
  183.      mp_ptr vp;
  184.      mp_size_t vsize;
  185.      mp_ptr up;
  186.      mp_size_t usize;
  187. #endif
  188. {
  189.   mp_ptr orig_vp = vp;
  190.   mp_size_t orig_vsize = vsize;
  191.   int binary_gcd_ctr;        /* Number of times binary gcd will execute.  */
  192.   TMP_DECL (marker);
  193.  
  194.   TMP_MARK (marker);
  195.  
  196.   /* Use accelerated algorithm if vsize is over ACCEL_THRESHOLD.
  197.      Two EXTRA limbs for U and V are required for kary reduction.  */
  198.   if (vsize > ACCEL_THRESHOLD)
  199.     {
  200.       unsigned long int vbitsize, d;
  201.       mp_ptr orig_up = up;
  202.       mp_size_t orig_usize = usize;
  203.       mp_ptr anchor_up = (mp_ptr) TMP_ALLOC ((usize + 2) * BYTES_PER_MP_LIMB);
  204.  
  205.       MPN_COPY (anchor_up, orig_up, usize);
  206.       up = anchor_up;
  207.  
  208.       count_leading_zeros (d, up[usize-1]);
  209.       d = usize * BITS_PER_MP_LIMB - d;
  210.       count_leading_zeros (vbitsize, vp[vsize-1]);
  211.       vbitsize = vsize * BITS_PER_MP_LIMB - vbitsize;
  212.       d = d - vbitsize + 1;
  213.  
  214.       /* Use bmod reduction to quickly discover whether V divides U.  */
  215.       up[usize++] = 0;                /* Insert leading zero.  */
  216.       mpn_bdivmod (up, up, usize, vp, vsize, d);
  217.  
  218.       /* Now skip U/V mod 2^d and any low zero limbs.  */
  219.       d /= BITS_PER_MP_LIMB, up += d, usize -= d;
  220.       while (usize != 0 && up[0] == 0)
  221.     up++, usize--;
  222.  
  223.       if (usize == 0)                /* GCD == ORIG_V.  */
  224.     goto done;
  225.  
  226.       vp = (mp_ptr) TMP_ALLOC ((vsize + 2) * BYTES_PER_MP_LIMB);
  227.       MPN_COPY (vp, orig_vp, vsize);
  228.  
  229.       do                    /* Main loop.  */
  230.     {
  231.       if (up[usize-1] & SIGN_BIT)        /* U < 0; take twos' compl. */
  232.         {
  233.           mp_size_t i;
  234.           anchor_up[0] = -up[0];
  235.           for (i = 1; i < usize; i++)
  236.         anchor_up[i] = ~up[i];
  237.           up = anchor_up;
  238.         }
  239.  
  240.       MPN_NORMALIZE_NOT_ZERO (up, usize);
  241.  
  242.       if ((up[0] & 1) == 0)            /* Result even; remove twos. */
  243.         {
  244.           unsigned long int r;
  245.           count_trailing_zeros (r, up[0]);
  246.           mpn_rshift (anchor_up, up, usize, r);
  247.           usize -= (anchor_up[usize-1] == 0);
  248.         }
  249.       else if (anchor_up != up)
  250.         MPN_COPY (anchor_up, up, usize);
  251.  
  252.       SWAP_MPN (anchor_up, usize, vp, vsize);
  253.       up = anchor_up;
  254.  
  255.       if (vsize <= 2)        /* Kary can't handle < 2 limbs and  */
  256.         break;            /* isn't efficient for == 2 limbs.  */
  257.  
  258.       d = vbitsize;
  259.       count_leading_zeros (vbitsize, vp[vsize-1]);
  260.       vbitsize = vsize * BITS_PER_MP_LIMB - vbitsize;
  261.       d = d - vbitsize + 1;
  262.  
  263.       if (d > BMOD_THRESHOLD)    /* Bmod reduction.  */
  264.         {
  265.           up[usize++] = 0;
  266.           mpn_bdivmod (up, up, usize, vp, vsize, d);
  267.           d /= BITS_PER_MP_LIMB, up += d, usize -= d;
  268.         }
  269.       else                /* Kary reduction.  */
  270.         {
  271.           mp_limb_t bp[2], cp[2];
  272.  
  273.           /* C <-- V/U mod 2^(2*BITS_PER_MP_LIMB).  */
  274.           cp[0] = vp[0], cp[1] = vp[1];
  275.           mpn_bdivmod (cp, cp, 2, up, 2, 2*BITS_PER_MP_LIMB);
  276.  
  277.           /* U <-- find_a (C)  *  U.  */
  278.           up[usize] = mpn_mul_1 (up, up, usize, find_a (cp));
  279.           usize++;
  280.  
  281.           /* B <-- A/C == U/V mod 2^(BITS_PER_MP_LIMB + 1).
  282.           bp[0] <-- U/V mod 2^BITS_PER_MP_LIMB and
  283.           bp[1] <-- ( (U - bp[0] * V)/2^BITS_PER_MP_LIMB ) / V mod 2 */
  284.           bp[0] = up[0], bp[1] = up[1];
  285.           mpn_bdivmod (bp, bp, 2, vp, 2, BITS_PER_MP_LIMB);
  286.           bp[1] &= 1;    /* Since V is odd, division is unnecessary.  */
  287.  
  288.           up[usize++] = 0;
  289.           if (bp[1])    /* B < 0: U <-- U + (-B)  * V.  */
  290.         {
  291.            mp_limb_t c = mpn_addmul_1 (up, vp, vsize, -bp[0]);
  292.            mpn_add_1 (up + vsize, up + vsize, usize - vsize, c);
  293.         }
  294.           else        /* B >= 0:  U <-- U - B * V.  */
  295.         {
  296.           mp_limb_t b = mpn_submul_1 (up, vp, vsize, bp[0]);
  297.           mpn_sub_1 (up + vsize, up + vsize, usize - vsize, b);
  298.         }
  299.  
  300.           up += 2, usize -= 2;  /* At least two low limbs are zero.  */
  301.         }
  302.  
  303.       /* Must remove low zero limbs before complementing.  */
  304.       while (usize != 0 && up[0] == 0)
  305.         up++, usize--;
  306.     }
  307.       while (usize);
  308.  
  309.       /* Compute GCD (ORIG_V, GCD (ORIG_U, V)).  Binary will execute twice.  */
  310.       up = orig_up, usize = orig_usize;
  311.       binary_gcd_ctr = 2;
  312.     }
  313.   else
  314.     binary_gcd_ctr = 1;
  315.  
  316.   /* Finish up with the binary algorithm.  Executes once or twice.  */
  317.   for ( ; binary_gcd_ctr--; up = orig_vp, usize = orig_vsize)
  318.     {
  319.       if (usize > 2)        /* First make U close to V in size.  */
  320.     {
  321.       unsigned long int vbitsize, d;
  322.       count_leading_zeros (d, up[usize-1]);
  323.       d = usize * BITS_PER_MP_LIMB - d;
  324.       count_leading_zeros (vbitsize, vp[vsize-1]);
  325.       vbitsize = vsize * BITS_PER_MP_LIMB - vbitsize;
  326.       d = d - vbitsize - 1;
  327.       if (d != -(unsigned long int)1 && d > 2)
  328.         {
  329.           mpn_bdivmod (up, up, usize, vp, vsize, d);  /* Result > 0.  */
  330.           d /= (unsigned long int)BITS_PER_MP_LIMB, up += d, usize -= d;
  331.         }
  332.     }
  333.  
  334.       /* Start binary GCD.  */
  335.       do
  336.     {
  337.       mp_size_t zeros;
  338.  
  339.       /* Make sure U is odd.  */
  340.       MPN_NORMALIZE (up, usize);
  341.       while (up[0] == 0)
  342.         up += 1, usize -= 1;
  343.       if ((up[0] & 1) == 0)
  344.         {
  345.           unsigned long int r;
  346.           count_trailing_zeros (r, up[0]);
  347.           mpn_rshift (up, up, usize, r);
  348.           usize -= (up[usize-1] == 0);
  349.         }
  350.  
  351.       /* Keep usize >= vsize.  */
  352.       if (usize < vsize)
  353.         SWAP_MPN (up, usize, vp, vsize);
  354.  
  355.       if (usize <= 2)                /* Double precision. */
  356.         {
  357.           if (vsize == 1)
  358.         vp[0] = mpn_gcd_1 (up, usize, vp[0]);
  359.           else
  360.         vsize = gcd_2 (vp, up);
  361.           break;                    /* Binary GCD done.  */
  362.         }
  363.  
  364.       /* Count number of low zero limbs of U - V.  */
  365.       for (zeros = 0; up[zeros] == vp[zeros] && ++zeros != vsize; )
  366.         continue;
  367.  
  368.       /* If U < V, swap U and V; in any case, subtract V from U.  */
  369.       if (zeros == vsize)                /* Subtract done.  */
  370.         up += zeros, usize -= zeros;
  371.       else if (usize == vsize)
  372.         {
  373.           mp_size_t size = vsize;
  374.           do
  375.         size--;
  376.           while (up[size] == vp[size]);
  377.           if (up[size] < vp[size])            /* usize == vsize.  */
  378.         SWAP_PTR (up, vp);
  379.           up += zeros, usize = size + 1 - zeros;
  380.           mpn_sub_n (up, up, vp + zeros, usize);
  381.         }
  382.       else
  383.         {
  384.           mp_size_t size = vsize - zeros;
  385.           up += zeros, usize -= zeros;
  386.           if (mpn_sub_n (up, up, vp + zeros, size))
  387.         {
  388.           while (up[size] == 0)            /* Propagate borrow. */
  389.             up[size++] = -(mp_limb_t)1;
  390.           up[size] -= 1;
  391.         }
  392.         }
  393.     }
  394.       while (usize);                    /* End binary GCD.  */
  395.     }
  396.  
  397. done:
  398.   if (vp != gp)
  399.     MPN_COPY (gp, vp, vsize);
  400.   TMP_FREE (marker);
  401.   return vsize;
  402. }
  403.