diff options
author | Kevin Ryde <user42@zip.com.au> | 2001-03-14 23:27:40 +0100 |
---|---|---|
committer | Kevin Ryde <user42@zip.com.au> | 2001-03-14 23:27:40 +0100 |
commit | 95d98dfc07bd04667064b40aa9da72f5605e1cd8 (patch) | |
tree | f242926f3d4f141c146d04937f05a9b59bd791c5 /mpz/cong.c | |
parent | ed34f9c4d0b64dc3081b98a10ccd07ac9536f58f (diff) | |
download | gmp-95d98dfc07bd04667064b40aa9da72f5605e1cd8.tar.gz |
* mpz/cong.c, mpz/cong_2exp.c, mpz/cong_ui.c: New files.
Diffstat (limited to 'mpz/cong.c')
-rw-r--r-- | mpz/cong.c | 174 |
1 files changed, 174 insertions, 0 deletions
diff --git a/mpz/cong.c b/mpz/cong.c new file mode 100644 index 000000000..d5f2d80ae --- /dev/null +++ b/mpz/cong.c @@ -0,0 +1,174 @@ +/* mpz_congruent_p -- test congruence of two mpz's */ + +/* +Copyright 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your +option) any later version. + +The GNU MP Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. +*/ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + + +/* For big divisors this code is only very slightly better than the user + doing a combination of mpz_sub and mpz_tdiv_r, but it's quite convenient, + and perhaps in the future can be improved, in similar ways to + mpn_divisible_p perhaps. + + The csize==1 / dsize==1 special case makes mpz_congruent_p as good as + mpz_congruent_ui_p on relevant operands, though such a combination + probably doesn't occur often. + + Alternatives: + + If c<d then it'd work to just form a%d and compare a and c (either as + a==c or a+c==d depending on the signs), but the saving from avoiding the + abs(a-c) calculation would be small compared to the division. + + Similarly if both a<d and c<d then it would work to just compare a and c + (a==c or a+c==d), but this isn't considered a particularly important case + and so isn't done for the moment. + + Low zero limbs on d could be stripped and the corresponding limbs of a + and c tested and skipped, but doing so would introduce a borrow when a + and c differ in sign and have non-zero skipped limbs. It doesn't seem + worth the complications to do this, since low zero limbs on d should + occur only rarely. */ + +int +mpz_congruent_p (mpz_srcptr a, mpz_srcptr c, mpz_srcptr d) +{ + mp_size_t asize, csize, dsize, xor; + mp_srcptr ap, cp, dp; + mp_ptr xp; + mp_limb_t alow, clow, dlow, dmask, r; + int result; + TMP_DECL (marker); + + dsize = SIZ(d); + if (dsize == 0) + DIVIDE_BY_ZERO; + dsize = ABS(dsize); + dp = PTR(d); + + if (ABSIZ(a) < ABSIZ(c)) + MPZ_SRCPTR_SWAP (a, c); + + asize = SIZ(a); + csize = SIZ(c); + xor = (asize ^ csize); + + asize = ABS(asize); + ap = PTR(a); + + if (csize == 0) + return mpn_divisible_p (ap, asize, dp, dsize); + + csize = ABS(csize); + cp = PTR(c); + + alow = ap[0]; + clow = cp[0]; + dlow = dp[0]; + + /* Check a==c mod low zero bits of dlow. This might catch a few cases of + a!=c quickly, and it helps the csize==1 special cases below. */ + dmask = LOW_ZEROS_MASK (dlow); + alow = (xor >= 0 ? alow : -alow); + if (((alow-clow) & dmask) != 0) + return 0; + + if (csize == 1) + { + if (dsize == 1) + { + cong_1: + if (xor < 0) + NEG_MOD (clow, clow, dlow); + + if (BELOW_THRESHOLD (asize, MODEXACT_1_ODD_THRESHOLD)) + { + r = mpn_mod_1 (ap, asize, dlow); + if (c < d) + return r == clow; + else + return r == (clow % dlow); + } + + if ((dlow & 1) == 0) + { + /* Strip low zero bits to get odd d required by modexact. If + d==e*2^n then a==c mod d if and only if both a==c mod e and + a==c mod 2^n, the latter having been done above. */ + unsigned twos; + count_trailing_zeros (twos, dlow); + dlow >>= twos; + } + + r = mpn_modexact_1c_odd (ap, asize, dlow, clow); + return r == 0 || r == dlow; + } + + /* dlow==0 is avoided since we don't want to bother handling extra low + zero bits if dsecond is even (would involve borrow if a,c differ in + sign and alow,clow!=0). */ + if (dsize == 2 && dlow != 0) + { + mp_limb_t dsecond = dp[1]; + + if (dsecond <= dmask) + { + unsigned twos; + count_trailing_zeros (twos, dlow); + dlow = (dlow >> twos) | (dsecond << (BITS_PER_MP_LIMB-twos)); + + /* dlow will be odd here, so the test for it even under cong_1 + is unnecessary, but the rest of that code is wanted. */ + goto cong_1; + } + } + } + + TMP_MARK (marker); + xp = TMP_ALLOC_LIMBS (asize+1); + + /* calculate abs(a-c) */ + if (xor >= 0) + { + /* same signs, subtract */ + if (asize > csize || mpn_cmp (ap, cp, asize) >= 0) + ASSERT_NOCARRY (mpn_sub (xp, ap, asize, cp, csize)); + else + ASSERT_NOCARRY (mpn_sub_n (xp, cp, ap, asize)); + MPN_NORMALIZE (xp, asize); + } + else + { + /* different signs, add */ + mp_limb_t carry; + carry = mpn_add (xp, ap, asize, cp, csize); + xp[asize] = carry; + asize += (carry != 0); + } + + result = mpn_divisible_p (xp, asize, dp, dsize); + + TMP_FREE (marker); + return result; +} |