source/libs/gmp/gmp-src/mpn/generic/divis.c

   1 /* mpn_divisible_p -- mpn by mpn divisibility test
   2
   3    THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY.  THEY'RE ALMOST
   4    CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
   5    FUTURE GNU MP RELEASES.
   6
   7 Copyright 2001, 2002, 2005, 2009, 2014 Free Software Foundation, Inc.
   8
   9 This file is part of the GNU MP Library.
  10
  11 The GNU MP Library is free software; you can redistribute it and/or modify
  12 it under the terms of either:
  13
  14   * the GNU Lesser General Public License as published by the Free
  15     Software Foundation; either version 3 of the License, or (at your
  16     option) any later version.
  17
  18 or
  19
  20   * the GNU General Public License as published by the Free Software
  21     Foundation; either version 2 of the License, or (at your option) any
  22     later version.
  23
  24 or both in parallel, as here.
  25
  26 The GNU MP Library is distributed in the hope that it will be useful, but
  27 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  28 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  29 for more details.
  30
  31 You should have received copies of the GNU General Public License and the
  32 GNU Lesser General Public License along with the GNU MP Library.  If not,
  33 see https://www.gnu.org/licenses/.  */
  34
  35 #include "gmp.h"
  36 #include "gmp-impl.h"
  37 #include "longlong.h"
  38
  39
  40 /* Determine whether A={ap,an} is divisible by D={dp,dn}.  Must have both
  41    operands normalized, meaning high limbs non-zero, except that an==0 is
  42    allowed.
  43
  44    There usually won't be many low zero bits on D, but the checks for this
  45    are fast and might pick up a few operand combinations, in particular they
  46    might reduce D to fit the single-limb mod_1/modexact_1 code.
  47
  48    Future:
  49
  50    Getting the remainder limb by limb would make an early exit possible on
  51    finding a non-zero.  This would probably have to be bdivmod style so
  52    there's no addback, but it would need a multi-precision inverse and so
  53    might be slower than the plain method (on small sizes at least).
  54
  55    When D must be normalized (shifted to low bit set), it's possible to
  56    suppress the bit-shifting of A down, as long as it's already been checked
  57    that A has at least as many trailing zero bits as D.  */
  58
  59 int
  60 mpn_divisible_p (mp_srcptr ap, mp_size_t an,
  61                  mp_srcptr dp, mp_size_t dn)
  62 {
  63   mp_limb_t  alow, dlow, dmask;
  64   mp_ptr     qp, rp, tp;
  65   mp_size_t  i;
  66   mp_limb_t di;
  67   unsigned  twos;
  68   TMP_DECL;
  69
  70   ASSERT (an >= 0);
  71   ASSERT (an == 0 || ap[an-1] != 0);
  72   ASSERT (dn >= 1);
  73   ASSERT (dp[dn-1] != 0);
  74   ASSERT_MPN (ap, an);
  75   ASSERT_MPN (dp, dn);
  76
  77   /* When a<d only a==0 is divisible.
  78      Notice this test covers all cases of an==0. */
  79   if (an < dn)
  80     return (an == 0);
  81
  82   /* Strip low zero limbs from d, requiring a==0 on those. */
  83   for (;;)
  84     {
  85       alow = *ap;
  86       dlow = *dp;
  87
  88       if (dlow != 0)
  89         break;
  90
  91       if (alow != 0)
  92         return 0;  /* a has fewer low zero limbs than d, so not divisible */
  93
  94       /* a!=0 and d!=0 so won't get to n==0 */
  95       an--; ASSERT (an >= 1);
  96       dn--; ASSERT (dn >= 1);
  97       ap++;
  98       dp++;
  99     }
 100
 101   /* a must have at least as many low zero bits as d */
 102   dmask = LOW_ZEROS_MASK (dlow);
 103   if ((alow & dmask) != 0)
 104     return 0;
 105
 106   if (dn == 1)
 107     {
 108       if (ABOVE_THRESHOLD (an, BMOD_1_TO_MOD_1_THRESHOLD))
 109         return mpn_mod_1 (ap, an, dlow) == 0;
 110
 111       count_trailing_zeros (twos, dlow);
 112       dlow >>= twos;
 113       return mpn_modexact_1_odd (ap, an, dlow) == 0;
 114     }
 115
 116   count_trailing_zeros (twos, dlow);
 117   if (dn == 2)
 118     {
 119       mp_limb_t  dsecond = dp[1];
 120       if (dsecond <= dmask)
 121         {
 122           dlow = (dlow >> twos) | (dsecond << (GMP_NUMB_BITS-twos));
 123           ASSERT_LIMB (dlow);
 124           return MPN_MOD_OR_MODEXACT_1_ODD (ap, an, dlow) == 0;
 125         }
 126     }
 127
 128   /* Should we compute Q = A * D^(-1) mod B^k,
 129                        R = A - Q * D  mod B^k
 130      here, for some small values of k?  Then check if R = 0 (mod B^k).  */
 131
 132   /* We could also compute A' = A mod T and D' = D mod P, for some
 133      P = 3 * 5 * 7 * 11 ..., and then check if any prime factor from P
 134      dividing D' also divides A'.  */
 135
 136   TMP_MARK;
 137
 138   TMP_ALLOC_LIMBS_2 (rp, an + 1,
 139                      qp, an - dn + 1); /* FIXME: Could we avoid this? */
 140
 141   if (twos != 0)
 142     {
 143       tp = TMP_ALLOC_LIMBS (dn);
 144       ASSERT_NOCARRY (mpn_rshift (tp, dp, dn, twos));
 145       dp = tp;
 146
 147       ASSERT_NOCARRY (mpn_rshift (rp, ap, an, twos));
 148     }
 149   else
 150     {
 151       MPN_COPY (rp, ap, an);
 152     }
 153   if (rp[an - 1] >= dp[dn - 1])
 154     {
 155       rp[an] = 0;
 156       an++;
 157     }
 158   else if (an == dn)
 159     {
 160       TMP_FREE;
 161       return 0;
 162     }
 163
 164   ASSERT (an > dn);             /* requirement of functions below */
 165
 166   if (BELOW_THRESHOLD (dn, DC_BDIV_QR_THRESHOLD) ||
 167       BELOW_THRESHOLD (an - dn, DC_BDIV_QR_THRESHOLD))
 168     {
 169       binvert_limb (di, dp[0]);
 170       mpn_sbpi1_bdiv_qr (qp, rp, an, dp, dn, -di);
 171       rp += an - dn;
 172     }
 173   else if (BELOW_THRESHOLD (dn, MU_BDIV_QR_THRESHOLD))
 174     {
 175       binvert_limb (di, dp[0]);
 176       mpn_dcpi1_bdiv_qr (qp, rp, an, dp, dn, -di);
 177       rp += an - dn;
 178     }
 179   else
 180     {
 181       tp = TMP_ALLOC_LIMBS (mpn_mu_bdiv_qr_itch (an, dn));
 182       mpn_mu_bdiv_qr (qp, rp, rp, an, dp, dn, tp);
 183     }
 184
 185   /* test for {rp,dn} zero or non-zero */
 186   i = 0;
 187   do
 188     {
 189       if (rp[i] != 0)
 190         {
 191           TMP_FREE;
 192           return 0;
 193         }
 194     }
 195   while (++i < dn);
 196
 197   TMP_FREE;
 198   return 1;
 199 }