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