contrib/gmp/mpz/cong.c

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