1/* mpz_congruent_p -- test congruence of two mpz's. 2 3Copyright 2001, 2002, 2005 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 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 14or 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 20or both in parallel, as here. 21 22The GNU MP Library is distributed in the hope that it will be useful, but 23WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 24or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 25for more details. 26 27You should have received copies of the GNU General Public License and the 28GNU Lesser General Public License along with the GNU MP Library. If not, 29see 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 60int 61mpz_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} 183