Annotation of OpenXM_contrib/gmp/mpn/generic/divis.c, Revision 1.1.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