contrib/gmp/mpz/cong_2exp.c

   1 /* mpz_congruent_2exp_p -- test congruence of mpz mod 2^n.
   2
   3 Copyright 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 int
  25 mpz_congruent_2exp_p (mpz_srcptr a, mpz_srcptr c, unsigned long d)
  26 {
  27   unsigned long  i, dlimbs, dbits;
  28   mp_ptr         ap, cp;
  29   mp_limb_t      dmask, alimb, climb, sum;
  30   mp_size_t      asize_signed, csize_signed, asize, csize;
  31
  32   if (ABSIZ(a) < ABSIZ(c))
  33     MPZ_SRCPTR_SWAP (a, c);
  34
  35   dlimbs = d / GMP_NUMB_BITS;
  36   dbits = d % GMP_NUMB_BITS;
  37   dmask = (CNST_LIMB(1) << dbits) - 1;
  38
  39   ap = PTR(a);
  40   cp = PTR(c);
  41
  42   asize_signed = SIZ(a);
  43   asize = ABS(asize_signed);
  44
  45   csize_signed = SIZ(c);
  46   csize = ABS(csize_signed);
  47
  48   if (csize_signed == 0)
  49     goto a_zeros;
  50
  51   if ((asize_signed ^ csize_signed) >= 0)
  52     {
  53       /* same signs, direct comparison */
  54
  55       /* a==c for limbs in common */
  56       if (mpn_cmp (ap, cp, MIN (csize, dlimbs)) != 0)
  57         return 0;
  58
  59       /* if that's all of dlimbs, then a==c for remaining bits */
  60       if (csize > dlimbs)
  61         return ((ap[dlimbs]-cp[dlimbs]) & dmask) == 0;
  62
  63     a_zeros:
  64       /* a remains, need all zero bits */
  65
  66       /* if d covers all of a and c, then must be exactly equal */
  67       if (asize <= dlimbs)
  68         return asize == csize;
  69
  70       /* whole limbs zero */
  71       for (i = csize; i < dlimbs; i++)
  72         if (ap[i] != 0)
  73           return 0;
  74
  75       /* partial limb zero */
  76       return (ap[dlimbs] & dmask) == 0;
  77     }
  78   else
  79     {
  80       /* different signs, negated comparison */
  81
  82       /* common low zero limbs, stopping at first non-zeros, which must
  83          match twos complement */
  84       i = 0;
  85       for (;;)
  86         {
  87           ASSERT (i < csize);  /* always have a non-zero limb on c */
  88           alimb = ap[i];
  89           climb = cp[i];
  90           sum = (alimb + climb) & GMP_NUMB_MASK;
  91
  92           if (i >= dlimbs)
  93             return (sum & dmask) == 0;
  94           i++;
  95
  96           /* require both zero, or first non-zeros as twos-complements */
  97           if (sum != 0)
  98             return 0;
  99
 100           if (alimb != 0)
 101             break;
 102         }
 103
 104       /* further limbs matching as ones-complement */
 105       for (;;)
 106         {
 107           if (i >= csize)
 108             break;
 109
 110           alimb = ap[i];
 111           climb = cp[i];
 112           sum = (alimb + climb + 1) & GMP_NUMB_MASK;
 113
 114           if (i >= dlimbs)
 115             return (sum & dmask) == 0;
 116
 117           if (sum != 0)
 118             return 0;
 119
 120           i++;
 121         }
 122
 123       /* no more c, so require all 1 bits in a */
 124
 125       if (asize < dlimbs)
 126         return 0;   /* not enough a */
 127
 128       /* whole limbs */
 129       for ( ; i < dlimbs; i++)
 130         if (ap[i] != GMP_NUMB_MAX)
 131           return 0;
 132
 133       /* if only whole limbs, no further fetches from a */
 134       if (dbits == 0)
 135         return 1;
 136
 137       /* need enough a */
 138       if (asize == dlimbs)
 139         return 0;
 140
 141       return ((ap[dlimbs]+1) & dmask) == 0;
 142     }
 143 }