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