Annotation of OpenXM_contrib/gmp/mpn/generic/divis.c, Revision 1.1
1.1 ! ohara 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 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 2.1 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; see the file COPYING.LIB. If not, write to
! 23: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 24: MA 02111-1307, USA. */
! 25:
! 26: #include "gmp.h"
! 27: #include "gmp-impl.h"
! 28: #include "longlong.h"
! 29:
! 30:
! 31: /* Determine whether {ap,asize} is divisible by {dp,dsize}. Must have both
! 32: operands normalized, meaning high limbs non-zero, except that asize==0 is
! 33: allowed.
! 34:
! 35: There usually won't be many low zero bits on d, but the checks for this
! 36: are fast and might pick up a few operand combinations, in particular they
! 37: might reduce d to fit the single-limb mod_1/modexact_1 code.
! 38:
! 39: Future:
! 40:
! 41: This is currently not much faster than the user doing an mpz_tdiv_r
! 42: and testing for a zero remainder, but hopefully it can be improved.
! 43:
! 44: mpn_bdivmod is one possibility, but it only trades udiv_qrnnd's for
! 45: multiplies, it won't save crossproducts the way it can in mpz_divexact.
! 46: Definitely worthwhile on small operands for most processors, but a
! 47: sub-quadratic version will be wanted before it can be used on all sizes.
! 48:
! 49: Getting the remainder limb by limb would make an early exit possible on
! 50: finding a non-zero. This would probably have to be bdivmod style so
! 51: there's no addback, but it would need a multi-precision inverse and so
! 52: might be slower than the plain method (on small sizes at least).
! 53:
! 54: When d must be normalized (shifted to high bit set), it's possible to
! 55: just append a low zero limb to "a" rather than bit-shifting as
! 56: mpn_tdiv_qr does internally, so long as it's already been checked that a
! 57: has at least as many trailing zeros bits as d. Or equivalently, pass
! 58: qxn==1 to mpn_tdiv_qr, if/when it accepts that.
! 59:
! 60: When called from mpz_congruent_p, {ap,asize} is a temporary which can be
! 61: destroyed. Maybe it'd be possible to get into mpn_tdiv_qr at a lower
! 62: level to save copying it, or maybe that function could accept rp==ap.
! 63:
! 64: Could use __attribute__ ((regparm (2))) on i386, so the parameters
! 65: wouldn't need extra stack when called from mpz_divisible_p, but a
! 66: pre-release gcc 3 didn't generate particularly good register juggling in
! 67: that case, so this isn't done for now. */
! 68:
! 69: int
! 70: mpn_divisible_p (mp_srcptr ap, mp_size_t asize,
! 71: mp_srcptr dp, mp_size_t dsize)
! 72: {
! 73: mp_limb_t alow, dlow, dmask;
! 74: mp_ptr qp, rp;
! 75: mp_size_t i;
! 76: TMP_DECL (marker);
! 77:
! 78: ASSERT (asize >= 0);
! 79: ASSERT (asize == 0 || ap[asize-1] != 0);
! 80: ASSERT (dsize >= 1);
! 81: ASSERT (dp[dsize-1] != 0);
! 82: ASSERT_MPN (ap, asize);
! 83: ASSERT_MPN (dp, dsize);
! 84:
! 85: /* When a<d only a==0 is divisible.
! 86: Notice this test covers all cases of asize==0. */
! 87: if (asize < dsize)
! 88: return (asize == 0);
! 89:
! 90: /* Strip low zero limbs from d, requiring a==0 on those. */
! 91: for (;;)
! 92: {
! 93: alow = *ap;
! 94: dlow = *dp;
! 95:
! 96: if (dlow != 0)
! 97: break;
! 98:
! 99: if (alow != 0)
! 100: return 0; /* a has fewer low zero limbs than d, so not divisible */
! 101:
! 102: /* a!=0 and d!=0 so won't get to size==0 */
! 103: asize--; ASSERT (asize >= 1);
! 104: dsize--; ASSERT (dsize >= 1);
! 105: ap++;
! 106: dp++;
! 107: }
! 108:
! 109: /* a must have at least as many low zero bits as d */
! 110: dmask = LOW_ZEROS_MASK (dlow);
! 111: if ((alow & dmask) != 0)
! 112: return 0;
! 113:
! 114: if (dsize == 1)
! 115: {
! 116: if (BELOW_THRESHOLD (asize, MODEXACT_1_ODD_THRESHOLD))
! 117: return mpn_mod_1 (ap, asize, dlow) == 0;
! 118:
! 119: if ((dlow & 1) == 0)
! 120: {
! 121: unsigned twos;
! 122: count_trailing_zeros (twos, dlow);
! 123: dlow >>= twos;
! 124: }
! 125: return mpn_modexact_1_odd (ap, asize, dlow) == 0;
! 126: }
! 127:
! 128: if (dsize == 2)
! 129: {
! 130: mp_limb_t dsecond = dp[1];
! 131: if (dsecond <= dmask)
! 132: {
! 133: unsigned twos;
! 134: count_trailing_zeros (twos, dlow);
! 135: dlow = (dlow >> twos) | (dsecond << (GMP_NUMB_BITS-twos));
! 136: ASSERT_LIMB (dlow);
! 137: return MPN_MOD_OR_MODEXACT_1_ODD (ap, asize, dlow) == 0;
! 138: }
! 139: }
! 140:
! 141: TMP_MARK (marker);
! 142:
! 143: rp = TMP_ALLOC_LIMBS (asize+1);
! 144: qp = rp + dsize;
! 145:
! 146: mpn_tdiv_qr (qp, rp, 0, ap, asize, dp, dsize);
! 147:
! 148: /* test for {rp,dsize} zero or non-zero */
! 149: i = 0;
! 150: do
! 151: {
! 152: if (rp[i] != 0)
! 153: {
! 154: TMP_FREE (marker);
! 155: return 0;
! 156: }
! 157: }
! 158: while (++i < dsize);
! 159:
! 160: TMP_FREE (marker);
! 161: return 1;
! 162: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to divis.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpn / generic