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