[BACK]Return to cong.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / gmp / mpz

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>