Annotation of OpenXM_contrib/gmp/mpz/cong_2exp.c, Revision 1.1.1.1
1.1 ohara 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 2.1 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; see the file COPYING.LIB. If not, write to
19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
20: MA 02111-1307, USA. */
21:
22: #include "gmp.h"
23: #include "gmp-impl.h"
24:
25:
26: int
27: mpz_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: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>