source: trunk/third/gmp/mpz/cong_2exp.c @ 18191

Revision 18191, 3.6 KB checked in by ghudson, 22 years ago (diff)
This commit was generated by cvs2svn to compensate for changes in r18190, which included commits to RCS files with non-trunk default branches.
Line 
1/* mpz_congruent_2exp_p -- test congruence of mpz mod 2^n.
2
3Copyright 2001, 2002 Free Software Foundation, Inc.
4
5This file is part of the GNU MP Library.
6
7The GNU MP Library is free software; you can redistribute it and/or modify
8it under the terms of the GNU Lesser General Public License as published by
9the Free Software Foundation; either version 2.1 of the License, or (at your
10option) any later version.
11
12The GNU MP Library is distributed in the hope that it will be useful, but
13WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
15License for more details.
16
17You should have received a copy of the GNU Lesser General Public License
18along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
19the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
20MA 02111-1307, USA. */
21
22#include "gmp.h"
23#include "gmp-impl.h"
24
25
26int
27mpz_congruent_2exp_p (mpz_srcptr a, mpz_srcptr c, unsigned long d)
28{
29  unsigned long  i, dlimbs, dbits;
30  mp_ptr         ap, cp;
31  mp_limb_t      dmask, alimb, climb, sum;
32  mp_size_t      asize_signed, csize_signed, asize, csize;
33
34  if (ABSIZ(a) < ABSIZ(c))
35    MPZ_SRCPTR_SWAP (a, c);
36
37  dlimbs = d / GMP_NUMB_BITS;
38  dbits = d % GMP_NUMB_BITS;
39  dmask = (CNST_LIMB(1) << dbits) - 1;
40
41  ap = PTR(a);
42  cp = PTR(c);
43
44  asize_signed = SIZ(a);
45  asize = ABS(asize_signed);
46
47  csize_signed = SIZ(c);
48  csize = ABS(csize_signed);
49
50  if (csize_signed == 0)
51    goto a_zeros;
52
53  if ((asize_signed ^ csize_signed) >= 0)
54    {
55      /* same signs, direct comparison */
56
57      /* a==c for limbs in common */
58      if (mpn_cmp (ap, cp, MIN (csize, dlimbs)) != 0)
59        return 0;
60
61      /* if that's all of dlimbs, then a==c for remaining bits */
62      if (csize > dlimbs)
63        return ((ap[dlimbs]-cp[dlimbs]) & dmask) == 0;
64
65    a_zeros:
66      /* a remains, need all zero bits */
67
68      /* if d covers all of a and c, then must be exactly equal */
69      if (asize <= dlimbs)
70        return asize == csize;
71
72      /* whole limbs zero */
73      for (i = csize; i < dlimbs; i++)
74        if (ap[i] != 0)
75          return 0;
76
77      /* partial limb zero */
78      return (ap[dlimbs] & dmask) == 0;
79    }
80  else
81    {
82      /* different signs, negated comparison */
83
84      /* common low zero limbs, stopping at first non-zeros, which must
85         match twos complement */
86      i = 0;
87      for (;;)
88        {
89          ASSERT (i < csize);  /* always have a non-zero limb on c */
90          alimb = ap[i];
91          climb = cp[i];
92          sum = (alimb + climb) & GMP_NUMB_MASK;
93
94          if (i >= dlimbs)
95            return (sum & dmask) == 0;
96          i++;
97
98          /* require both zero, or first non-zeros as twos-complements */
99          if (sum != 0)
100            return 0;
101
102          if (alimb != 0)
103            break;
104        }
105
106      /* further limbs matching as ones-complement */
107      for (;;)
108        {
109          if (i >= csize)
110            break;
111
112          alimb = ap[i];
113          climb = cp[i];
114          sum = (alimb + climb + 1) & GMP_NUMB_MASK;
115
116          if (i >= dlimbs)
117            return (sum & dmask) == 0;
118
119          if (sum != 0)
120            return 0;
121
122          i++;
123        }
124
125      /* no more c, so require all 1 bits in a */
126
127      if (asize < dlimbs)
128        return 0;   /* not enough a */
129
130      /* whole limbs */
131      for ( ; i < dlimbs; i++)
132        if (ap[i] != GMP_NUMB_MAX)
133          return 0;
134
135      /* if only whole limbs, no further fetches from a */
136      if (dbits == 0)
137        return 1;
138
139      /* need enough a */
140      if (asize == dlimbs)
141        return 0;
142
143      return ((ap[dlimbs]+1) & dmask) == 0;
144    }
145}
Note: See TracBrowser for help on using the repository browser.