* mpz/hamdist.c: Support neg/neg operands.

1 files changed, 128 insertions, 20 deletions

diff --git a/mpz/hamdist.c b/mpz/hamdist.c
index 06ca8d2e0..f8f12bf62 100644
--- a/mpz/hamdist.c
+++ b/mpz/hamdist.c
@@ -1,7 +1,4 @@
-/* mpz_hamdist(mpz_ptr op1, mpz_ptr op2) -- Compute the hamming distance
-   between OP1 and OP2.  If one of the operands is negative, return ~0.  (We
-   could make the function well-defined when both operands are negative, but
-   that would probably not be worth the trouble.
+/* mpz_hamdist -- calculate hamming distance.
 
 Copyright 1994, 1996, 2001 Free Software Foundation, Inc.
 
@@ -26,32 +23,143 @@ MA 02111-1307, USA. */
 #include "gmp-impl.h"
 
 
-unsigned long int
+unsigned long
 mpz_hamdist (mpz_srcptr u, mpz_srcptr v)
 {
-  mp_srcptr up, vp;
-  mp_size_t usize, vsize;
-  unsigned long int count;
+  mp_srcptr      up, vp;
+  mp_size_t      usize, vsize;
+  unsigned long  count;
 
   usize = SIZ(u);
   vsize = SIZ(v);
 
-  if ((usize | vsize) < 0)
-    return ~ (unsigned long int) 0;
-
   up = PTR(u);
   vp = PTR(v);
 
-  if (usize < vsize)
-    MPN_SRCPTR_SWAP (up,usize, vp,vsize);
+  if (usize >= 0)
+    {
+      if (vsize < 0)
+        return ~ (unsigned long) 0;
+
+      /* positive/positive */
+
+      if (usize < vsize)
+        MPN_SRCPTR_SWAP (up,usize, vp,vsize);
+
+      count = 0;
+      if (vsize != 0)
+        count = mpn_hamdist (up, vp, vsize);
+
+      usize -= vsize;
+      if (usize != 0)
+        count += mpn_popcount (up + vsize, usize);
+
+      return count;
+    }
+  else
+    {
+      mp_limb_t  ulimb, vlimb;
+      mp_size_t  old_vsize, step;
+
+      if (vsize >= 0)
+        return ~ (unsigned long) 0;
+
+      /* negative/negative */
+
+      usize = -usize;
+      vsize = -vsize;
+
+      /* skip common low zeros */
+      for (;;)
+        {
+          ASSERT (usize > 0);
+          ASSERT (vsize > 0);
+
+          usize--;
+          vsize--;
+
+          ulimb = *up++;
+          vlimb = *vp++;
+
+          if (ulimb != 0)
+            break;
+
+          if (vlimb != 0)
+            {
+              MPN_SRCPTR_SWAP (up,usize, vp,vsize);
+              ulimb = vlimb;
+              vlimb = 0;
+              break;
+            }
+        }
+
+      /* twos complement for non-zero limbs (ulimb is non-zero, but vlimb
+         might be zero) */
+      ulimb = -ulimb;
+      vlimb = -vlimb;
+      popc_limb (count, ulimb ^ vlimb);
+
+      if (vlimb == 0)
+        {
+          unsigned long  twoscount;
+
+          /* first non-zero of v */
+          old_vsize = vsize;
+          do
+            {
+              ASSERT (vsize > 0);
+              vsize--;
+              vlimb = *vp++;
+            }
+          while (vlimb == 0);
 
-  count = 0;
-  if (vsize != 0)
-    count = mpn_hamdist (up, vp, vsize);
+          /* part of u corresponding to skipped v zeros */
+          step = old_vsize - vsize - 1;
+          count += step * BITS_PER_MP_LIMB;
+          step = MIN (step, usize);
+          if (step != 0)
+            {
+              count -= mpn_popcount (up, step);
+              usize -= step;
+              up += step;
+            }
 
-  usize -= vsize;
-  if (usize != 0)
-    count += mpn_popcount (up + vsize, usize);
+          /* First non-zero vlimb as twos complement, xor with ones
+             complement ulimb.  Note -v^(~0^u) == (v-1)^u. */
+          vlimb--;
+          if (usize != 0)
+            {
+              usize--;
+              vlimb ^= *up++;
+            }
+          popc_limb (twoscount, vlimb);
+          count += twoscount;
+        }
+      
+      /* Overlapping part of u and v, if any.  Ones complement both, so just
+         plain hamdist. */
+      step = MIN (usize, vsize);
+      if (step != 0)
+        {
+          count += mpn_hamdist (up, vp, vsize);
+          usize -= step;
+          vsize -= step;
+          up += step;
+          vp += step;
+        }
 
-  return count;
+      /* Remaining high part of u or v, if any, ones complement. */
+      if (usize != 0)
+        {
+        remaining:
+          count += usize * BITS_PER_MP_LIMB - mpn_popcount (up, usize);
+        }
+      else if (vsize != 0)
+        {
+          up = vp;
+          usize = vsize;
+          goto remaining;
+        }
+      return count;
+    }
 }