OpenXM_contrib/gmp/mpn/generic/gcd_1.c - diff

Return to gcd_1.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpn / generic

Diff for /OpenXM_contrib/gmp/mpn/generic/Attic/gcd_1.c between version 1.1.1.2 and 1.1.1.3

-version 1.1.1.2, 2000/09/09 14:12:25
+version 1.1.1.3, 2003/08/25 16:06:20
 Line 1
 Line 1
 Line 1
- /* mpn_gcd_1 --
+ /* mpn_gcd_1 -- mpn and limb greatest common divisor.
- Copyright (C) 1994, 1996, 2000 Free Software Foundation, Inc.
+ Copyright 1994, 1996, 2000, 2001 Free Software Foundation, Inc.
  This file is part of the GNU MP Library.
 Line 23  MA 02111-1307, USA. */
 Line 23  MA 02111-1307, USA. */
 Line 23  MA 02111-1307, USA. */
  #include "gmp-impl.h"
  #include "longlong.h"
  /* Does not work for U == 0 or V == 0.  It would be tough to make it work for
-    V == 0 since gcd(x,0) = x, and U does not generally fit in an mp_limb_t.  */
+    V == 0 since gcd(x,0) = x, and U does not generally fit in an mp_limb_t.
+    The threshold for doing u%v when size==1 will vary by CPU according to
+    the speed of a division and the code generated for the main loop.  Any
+    tuning for this is left to a CPU specific implementation.  */
  mp_limb_t
- #if __STDC__
  mpn_gcd_1 (mp_srcptr up, mp_size_t size, mp_limb_t vlimb)
- #else
- mpn_gcd_1 (up, size, vlimb)
-      mp_srcptr up;
-      mp_size_t size;
-      mp_limb_t vlimb;
- #endif
  {
-   mp_limb_t ulimb;
+   mp_limb_t      ulimb;
-   unsigned long int u_low_zero_bits, v_low_zero_bits;
+   unsigned long  zero_bits, u_low_zero_bits;
+   ASSERT (size >= 1);
+   ASSERT (vlimb != 0);
+   ASSERT_MPN_NONZERO_P (up, size);
+   ulimb = up[0];
+   /* Need vlimb odd for modexact, want it odd to get common zeros. */
+   count_trailing_zeros (zero_bits, vlimb);
+   vlimb >>= zero_bits;
    if (size > 1)
      {
-       ulimb = mpn_mod_1 (up, size, vlimb);
+       /* Must get common zeros before the mod reduction.  If ulimb==0 then
+          vlimb already gives the common zeros.  */
+       if (ulimb != 0)
+         {
+           count_trailing_zeros (u_low_zero_bits, ulimb);
+           zero_bits = MIN (zero_bits, u_low_zero_bits);
+         }
+       ulimb = MPN_MOD_OR_MODEXACT_1_ODD (up, size, vlimb);
        if (ulimb == 0)
-         return vlimb;
+         goto done;
+       goto strip_u_maybe;
      }
-   else
-     ulimb = up[0];
-   /*  Need to eliminate low zero bits.  */
+   /* size==1, so up[0]!=0 */
    count_trailing_zeros (u_low_zero_bits, ulimb);
    ulimb >>= u_low_zero_bits;
+   zero_bits = MIN (zero_bits, u_low_zero_bits);
-   count_trailing_zeros (v_low_zero_bits, vlimb);
+   /* make u bigger */
-   vlimb >>= v_low_zero_bits;
+   if (vlimb > ulimb)
+     MP_LIMB_T_SWAP (ulimb, vlimb);
+   /* if u is much bigger than v, reduce using a division rather than
+      chipping away at it bit-by-bit */
+   if ((ulimb >> 16) > vlimb)
+     {
+       ulimb %= vlimb;
+       if (ulimb == 0)
+         goto done;
+       goto strip_u_maybe;
+     }
    while (ulimb != vlimb)
      {
+       ASSERT (ulimb & 1);
+       ASSERT (vlimb & 1);
        if (ulimb > vlimb)
          {
            ulimb -= vlimb;
            do
-             ulimb >>= 1;
+             {
+               ulimb >>= 1;
+               ASSERT (ulimb != 0);
+             strip_u_maybe:
+               ;
+             }
            while ((ulimb & 1) == 0);
          }
        else /*  vlimb > ulimb.  */
          {
            vlimb -= ulimb;
            do
-             vlimb >>= 1;
+             {
+               vlimb >>= 1;
+               ASSERT (vlimb != 0);
+             }
            while ((vlimb & 1) == 0);
          }
      }
-   return  ulimb << MIN (u_low_zero_bits, v_low_zero_bits);
+  done:
+   return vlimb << zero_bits;
  }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>