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