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