1/* mpz_hamdist -- calculate hamming distance.
2
3Copyright 1994, 1996, 2001, 2002, 2009, 2010, 2011 Free Software Foundation,
4Inc.
5
6This file is part of the GNU MP Library.
7
8The GNU MP Library is free software; you can redistribute it and/or modify
9it under the terms of the GNU Lesser General Public License as published by
10the Free Software Foundation; either version 3 of the License, or (at your
11option) any later version.
12
13The GNU MP Library is distributed in the hope that it will be useful, but
14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16License for more details.
17
18You should have received a copy of the GNU Lesser General Public License
19along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
20
21#include "gmp.h"
22#include "gmp-impl.h"
23
24
25mp_bitcnt_t
26mpz_hamdist (mpz_srcptr u, mpz_srcptr v) __GMP_NOTHROW
27{
28  mp_srcptr      up, vp;
29  mp_size_t      usize, vsize;
30  mp_bitcnt_t    count;
31
32  usize = SIZ(u);
33  vsize = SIZ(v);
34
35  up = PTR(u);
36  vp = PTR(v);
37
38  if (usize >= 0)
39    {
40      if (vsize < 0)
41	return ~ (mp_bitcnt_t) 0;
42
43      /* positive/positive */
44
45      if (usize < vsize)
46	MPN_SRCPTR_SWAP (up,usize, vp,vsize);
47
48      count = 0;
49      if (vsize != 0)
50	count = mpn_hamdist (up, vp, vsize);
51
52      usize -= vsize;
53      if (usize != 0)
54	count += mpn_popcount (up + vsize, usize);
55
56      return count;
57    }
58  else
59    {
60      mp_limb_t  ulimb, vlimb;
61      mp_size_t  old_vsize, step;
62
63      if (vsize >= 0)
64	return ~ (mp_bitcnt_t) 0;
65
66      /* negative/negative */
67
68      usize = -usize;
69      vsize = -vsize;
70
71      /* skip common low zeros */
72      for (;;)
73	{
74	  ASSERT (usize > 0);
75	  ASSERT (vsize > 0);
76
77	  usize--;
78	  vsize--;
79
80	  ulimb = *up++;
81	  vlimb = *vp++;
82
83	  if (ulimb != 0)
84	    break;
85
86	  if (vlimb != 0)
87	    {
88	      MPN_SRCPTR_SWAP (up,usize, vp,vsize);
89	      ulimb = vlimb;
90	      vlimb = 0;
91	      break;
92	    }
93	}
94
95      /* twos complement first non-zero limbs (ulimb is non-zero, but vlimb
96	 might be zero) */
97      ulimb = -ulimb;
98      vlimb = -vlimb;
99      popc_limb (count, (ulimb ^ vlimb) & GMP_NUMB_MASK);
100
101      if (vlimb == 0)
102	{
103	  mp_bitcnt_t  twoscount;
104
105	  /* first non-zero of v */
106	  old_vsize = vsize;
107	  do
108	    {
109	      ASSERT (vsize > 0);
110	      vsize--;
111	      vlimb = *vp++;
112	    }
113	  while (vlimb == 0);
114
115	  /* part of u corresponding to skipped v zeros */
116	  step = old_vsize - vsize - 1;
117	  count += step * GMP_NUMB_BITS;
118	  step = MIN (step, usize);
119	  if (step != 0)
120	    {
121	      count -= mpn_popcount (up, step);
122	      usize -= step;
123	      up += step;
124	    }
125
126	  /* First non-zero vlimb as twos complement, xor with ones
127	     complement ulimb.  Note -v^(~0^u) == (v-1)^u. */
128	  vlimb--;
129	  if (usize != 0)
130	    {
131	      usize--;
132	      vlimb ^= *up++;
133	    }
134	  popc_limb (twoscount, vlimb);
135	  count += twoscount;
136	}
137
138      /* Overlapping part of u and v, if any.  Ones complement both, so just
139	 plain hamdist. */
140      step = MIN (usize, vsize);
141      if (step != 0)
142	{
143	  count += mpn_hamdist (up, vp, step);
144	  usize -= step;
145	  vsize -= step;
146	  up += step;
147	  vp += step;
148	}
149
150      /* Remaining high part of u or v, if any, ones complement but xor
151	 against all ones in the other, so plain popcount. */
152      if (usize != 0)
153	{
154	remaining:
155	  count += mpn_popcount (up, usize);
156	}
157      else if (vsize != 0)
158	{
159	  up = vp;
160	  usize = vsize;
161	  goto remaining;
162	}
163      return count;
164    }
165}
166