contrib/gmp/mpz/hamdist.c

   1 /* mpz_hamdist -- calculate hamming distance.
   2
   3 Copyright 1994, 1996, 2001, 2002, 2009, 2010, 2011 Free Software Foundation,
   4 Inc.
   5
   6 This file is part of the GNU MP Library.
   7
   8 The GNU MP Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 3 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MP Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along 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
  25 mp_bitcnt_t
  26 mpz_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 }