Line | Count | Source (jump to first uncovered line) |
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 | 0 | { |
63 | 0 | mp_size_t asize, csize, dsize, sign; |
64 | 0 | mp_srcptr ap, cp, dp; |
65 | 0 | mp_ptr xp; |
66 | 0 | mp_limb_t alow, clow, dlow, dmask, r; |
67 | 0 | int result; |
68 | 0 | TMP_DECL; |
69 | |
|
70 | 0 | dsize = SIZ(d); |
71 | 0 | if (UNLIKELY (dsize == 0)) |
72 | 0 | return (mpz_cmp (a, c) == 0); |
73 | | |
74 | 0 | dsize = ABS(dsize); |
75 | 0 | dp = PTR(d); |
76 | |
|
77 | 0 | if (ABSIZ(a) < ABSIZ(c)) |
78 | 0 | MPZ_SRCPTR_SWAP (a, c); |
79 | |
|
80 | 0 | asize = SIZ(a); |
81 | 0 | csize = SIZ(c); |
82 | 0 | sign = (asize ^ csize); |
83 | |
|
84 | 0 | asize = ABS(asize); |
85 | 0 | ap = PTR(a); |
86 | |
|
87 | 0 | if (csize == 0) |
88 | 0 | return mpn_divisible_p (ap, asize, dp, dsize); |
89 | | |
90 | 0 | csize = ABS(csize); |
91 | 0 | cp = PTR(c); |
92 | |
|
93 | 0 | alow = ap[0]; |
94 | 0 | clow = cp[0]; |
95 | 0 | 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 | 0 | dmask = LOW_ZEROS_MASK (dlow) & GMP_NUMB_MASK; |
100 | 0 | alow = (sign >= 0 ? alow : -alow); |
101 | 0 | if (((alow-clow) & dmask) != 0) |
102 | 0 | return 0; |
103 | | |
104 | 0 | if (csize == 1) |
105 | 0 | { |
106 | 0 | if (dsize == 1) |
107 | 0 | { |
108 | 0 | cong_1: |
109 | 0 | if (sign < 0) |
110 | 0 | NEG_MOD (clow, clow, dlow); |
111 | |
|
112 | 0 | if (ABOVE_THRESHOLD (asize, BMOD_1_TO_MOD_1_THRESHOLD)) |
113 | 0 | { |
114 | 0 | r = mpn_mod_1 (ap, asize, dlow); |
115 | 0 | if (clow < dlow) |
116 | 0 | return r == clow; |
117 | 0 | else |
118 | 0 | return r == (clow % dlow); |
119 | 0 | } |
120 | | |
121 | 0 | if ((dlow & 1) == 0) |
122 | 0 | { |
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 | 0 | unsigned twos; |
127 | 0 | count_trailing_zeros (twos, dlow); |
128 | 0 | dlow >>= twos; |
129 | 0 | } |
130 | |
|
131 | 0 | r = mpn_modexact_1c_odd (ap, asize, dlow, clow); |
132 | 0 | return r == 0 || r == dlow; |
133 | 0 | } |
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 | 0 | if (dsize == 2 && dlow != 0) |
139 | 0 | { |
140 | 0 | mp_limb_t dsecond = dp[1]; |
141 | |
|
142 | 0 | if (dsecond <= dmask) |
143 | 0 | { |
144 | 0 | unsigned twos; |
145 | 0 | count_trailing_zeros (twos, dlow); |
146 | 0 | dlow = (dlow >> twos) | (dsecond << (GMP_NUMB_BITS-twos)); |
147 | 0 | 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 | 0 | goto cong_1; |
152 | 0 | } |
153 | 0 | } |
154 | 0 | } |
155 | | |
156 | 0 | TMP_MARK; |
157 | 0 | xp = TMP_ALLOC_LIMBS (asize+1); |
158 | | |
159 | | /* calculate abs(a-c) */ |
160 | 0 | if (sign >= 0) |
161 | 0 | { |
162 | | /* same signs, subtract */ |
163 | 0 | if (asize > csize || mpn_cmp (ap, cp, asize) >= 0) |
164 | 0 | ASSERT_NOCARRY (mpn_sub (xp, ap, asize, cp, csize)); |
165 | 0 | else |
166 | 0 | ASSERT_NOCARRY (mpn_sub_n (xp, cp, ap, asize)); |
167 | 0 | MPN_NORMALIZE (xp, asize); |
168 | 0 | } |
169 | 0 | else |
170 | 0 | { |
171 | | /* different signs, add */ |
172 | 0 | mp_limb_t carry; |
173 | 0 | carry = mpn_add (xp, ap, asize, cp, csize); |
174 | 0 | xp[asize] = carry; |
175 | 0 | asize += (carry != 0); |
176 | 0 | } |
177 | |
|
178 | 0 | result = mpn_divisible_p (xp, asize, dp, dsize); |
179 | |
|
180 | 0 | TMP_FREE; |
181 | 0 | return result; |
182 | 0 | } |