contrib/gmp/mpz/hamdist.c

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