Annotation of OpenXM_contrib/gmp/mpz/divegcd.c, Revision 1.1.1.1
1.1 ohara 1: /* mpz_divexact_gcd -- exact division optimized for GCDs.
2:
3: THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE AND ARE ALMOST CERTAIN TO
4: BE SUBJECT TO INCOMPATIBLE CHANGES IN FUTURE GNU MP RELEASES.
5:
6: Copyright 2000 Free Software Foundation, Inc.
7:
8: This file is part of the GNU MP Library.
9:
10: The GNU MP Library is free software; you can redistribute it and/or modify
11: it under the terms of the GNU Lesser General Public License as published by
12: the Free Software Foundation; either version 2.1 of the License, or (at your
13: option) any later version.
14:
15: The GNU MP Library is distributed in the hope that it will be useful, but
16: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
17: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18: License for more details.
19:
20: You should have received a copy of the GNU Lesser General Public License
21: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
22: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
23: MA 02111-1307, USA. */
24:
25: #include "gmp.h"
26: #include "gmp-impl.h"
27: #include "longlong.h"
28:
29:
30: /* Set q to a/d, expecting d to be from a GCD and therefore usually small.
31:
32: The distribution of GCDs of random numbers can be found in Knuth volume 2
33: section 4.5.2 theorem D.
34:
35: GCD chance
36: 1 60.7%
37: 2^k 20.2% (1<=k<32)
38: 3*2^k 9.0% (1<=k<32)
39: other 10.1%
40:
41: Only the low limb is examined for optimizations, since GCDs bigger than
42: 2^32 (or 2^64) will occur very infrequently.
43:
44: Future: This could change to an mpn_divexact_gcd, possibly partly
45: inlined, if/when the relevant mpq functions change to an mpn based
46: implementation. */
47:
48:
49: static void
50: mpz_divexact_by3 (mpz_ptr q, mpz_srcptr a)
51: {
52: mp_size_t size = SIZ(a);
53: if (size == 0)
54: {
55: SIZ(q) = 0;
56: return;
57: }
58: else
59: {
60: mp_size_t abs_size = ABS(size);
61: mp_ptr qp;
62:
63: MPZ_REALLOC (q, abs_size);
64:
65: qp = PTR(q);
66: mpn_divexact_by3 (qp, PTR(a), abs_size);
67:
68: abs_size -= (qp[abs_size-1] == 0);
69: SIZ(q) = (size>0 ? abs_size : -abs_size);
70: }
71: }
72:
73: void
74: mpz_divexact_gcd (mpz_ptr q, mpz_srcptr a, mpz_srcptr d)
75: {
76: ASSERT (mpz_sgn (d) > 0);
77:
78: if (SIZ(d) == 1)
79: {
80: mp_limb_t dl = PTR(d)[0];
81: int twos;
82:
83: if (dl == 1)
84: {
85: if (q != a)
86: mpz_set (q, a);
87: return;
88: }
89: if (dl == 3)
90: {
91: mpz_divexact_by3 (q, a);
92: return;
93: }
94:
95: count_trailing_zeros (twos, dl);
96: dl >>= twos;
97:
98: if (dl == 1)
99: {
100: mpz_tdiv_q_2exp (q, a, twos);
101: return;
102: }
103: if (dl == 3)
104: {
105: mpz_tdiv_q_2exp (q, a, twos);
106: mpz_divexact_by3 (q, q);
107: return;
108: }
109: }
110:
111: mpz_divexact (q, a, d);
112: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>