1 /* mpz_congruent_p -- test congruence of two mpz's.
2
3 Copyright 2001, 2002, 2005 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 #include "longlong.h"
33
34
35 /* For big divisors this code is only very slightly better than the user
36 doing a combination of mpz_sub and mpz_tdiv_r, but it's quite convenient,
37 and perhaps in the future can be improved, in similar ways to
38 mpn_divisible_p perhaps.
39
40 The csize==1 / dsize==1 special case makes mpz_congruent_p as good as
41 mpz_congruent_ui_p on relevant operands, though such a combination
42 probably doesn't occur often.
43
44 Alternatives:
45
46 If c<d then it'd work to just form a%d and compare a and c (either as
47 a==c or a+c==d depending on the signs), but the saving from avoiding the
48 abs(a-c) calculation would be small compared to the division.
49
50 Similarly if both a<d and c<d then it would work to just compare a and c
51 (a==c or a+c==d), but this isn't considered a particularly important case
52 and so isn't done for the moment.
53
54 Low zero limbs on d could be stripped and the corresponding limbs of a
55 and c tested and skipped, but doing so would introduce a borrow when a
56 and c differ in sign and have non-zero skipped limbs. It doesn't seem
57 worth the complications to do this, since low zero limbs on d should
58 occur only rarely. */
59
60 int
61 mpz_congruent_p (mpz_srcptr a, mpz_srcptr c, mpz_srcptr d)
62 {
63 mp_size_t asize, csize, dsize, sign;
64 mp_srcptr ap, cp, dp;
65 mp_ptr xp;
66 mp_limb_t alow, clow, dlow, dmask, r;
67 int result;
68 TMP_DECL;
69
70 dsize = SIZ(d);
71 if (UNLIKELY (dsize == 0))
72 return (mpz_cmp (a, c) == 0);
73
74 dsize = ABS(dsize);
75 dp = PTR(d);
76
77 if (ABSIZ(a) < ABSIZ(c))
78 MPZ_SRCPTR_SWAP (a, c);
79
80 asize = SIZ(a);
81 csize = SIZ(c);
82 sign = (asize ^ csize);
83
84 asize = ABS(asize);
85 ap = PTR(a);
86
87 if (csize == 0)
88 return mpn_divisible_p (ap, asize, dp, dsize);
89
90 csize = ABS(csize);
91 cp = PTR(c);
92
93 alow = ap[0];
94 clow = cp[0];
95 dlow = dp[0];
96
97 /* Check a==c mod low zero bits of dlow. This might catch a few cases of
98 a!=c quickly, and it helps the csize==1 special cases below. */
99 dmask = LOW_ZEROS_MASK (dlow) & GMP_NUMB_MASK;
100 alow = (sign >= 0 ? alow : -alow);
101 if (((alow-clow) & dmask) != 0)
102 return 0;
103
104 if (csize == 1)
105 {
106 if (dsize == 1)
107 {
108 cong_1:
109 if (sign < 0)
110 NEG_MOD (clow, clow, dlow);
111
112 if (ABOVE_THRESHOLD (asize, BMOD_1_TO_MOD_1_THRESHOLD))
113 {
114 r = mpn_mod_1 (ap, asize, dlow);
115 if (clow < dlow)
116 return r == clow;
117 else
118 return r == (clow % dlow);
119 }
120
121 if ((dlow & 1) == 0)
122 {
123 /* Strip low zero bits to get odd d required by modexact. If
124 d==e*2^n then a==c mod d if and only if both a==c mod e and
125 a==c mod 2^n, the latter having been done above. */
126 unsigned twos;
127 count_trailing_zeros (twos, dlow);
128 dlow >>= twos;
129 }
130
131 r = mpn_modexact_1c_odd (ap, asize, dlow, clow);
132 return r == 0 || r == dlow;
133 }
134
135 /* dlow==0 is avoided since we don't want to bother handling extra low
136 zero bits if dsecond is even (would involve borrow if a,c differ in
137 sign and alow,clow!=0). */
138 if (dsize == 2 && dlow != 0)
139 {
140 mp_limb_t dsecond = dp[1];
141
142 if (dsecond <= dmask)
143 {
144 unsigned twos;
145 count_trailing_zeros (twos, dlow);
146 dlow = (dlow >> twos) | (dsecond << (GMP_NUMB_BITS-twos));
147 ASSERT_LIMB (dlow);
148
149 /* dlow will be odd here, so the test for it even under cong_1
150 is unnecessary, but the rest of that code is wanted. */
151 goto cong_1;
152 }
153 }
154 }
155
156 TMP_MARK;
157 xp = TMP_ALLOC_LIMBS (asize+1);
158
159 /* calculate abs(a-c) */
160 if (sign >= 0)
161 {
162 /* same signs, subtract */
163 if (asize > csize || mpn_cmp (ap, cp, asize) >= 0)
164 ASSERT_NOCARRY (mpn_sub (xp, ap, asize, cp, csize));
165 else
166 ASSERT_NOCARRY (mpn_sub_n (xp, cp, ap, asize));
167 MPN_NORMALIZE (xp, asize);
168 }
169 else
170 {
171 /* different signs, add */
172 mp_limb_t carry;
173 carry = mpn_add (xp, ap, asize, cp, csize);
174 xp[asize] = carry;
175 asize += (carry != 0);
176 }
177
178 result = mpn_divisible_p (xp, asize, dp, dsize);
179
180 TMP_FREE;
181 return result;
182 }