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