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, 2005, 2011, 2012 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 either:
12
13 * the GNU Lesser General Public License as published by the Free
14 Software Foundation; either version 3 of the License, or (at your
15 option) any later version.
16
17 or
18
19 * the GNU General Public License as published by the Free Software
20 Foundation; either version 2 of the License, or (at your option) any
21 later version.
22
23 or both in parallel, as here.
24
25 The GNU MP Library is distributed in the hope that it will be useful, but
26 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
27 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
28 for more details.
29
30 You should have received copies of the GNU General Public License and the
31 GNU Lesser General Public License along with the GNU MP Library. If not,
32 see https://www.gnu.org/licenses/. */
33
34 #include "gmp-impl.h"
35 #include "longlong.h"
36
37
38 /* Set q to a/d, expecting d to be from a GCD and therefore usually small.
39
40 The distribution of GCDs of random numbers can be found in Knuth volume 2
41 section 4.5.2 theorem D.
42
43 GCD chance
44 1 60.8%
45 2^k 20.2% (1<=k<32)
46 3*2^k 9.0% (1<=k<32)
47 other 10.1%
48
49 Only the low limb is examined for optimizations, since GCDs bigger than
50 2^32 (or 2^64) will occur very infrequently.
51
52 Future: This could change to an mpn_divexact_gcd, possibly partly
53 inlined, if/when the relevant mpq functions change to an mpn based
54 implementation. */
55
56
57 #if GMP_NUMB_BITS % 2 == 0
58 static void
59 mpz_divexact_by3 (mpz_ptr q, mpz_srcptr a)
60 {
61 mp_size_t size = SIZ(a);
62 mp_size_t abs_size = ABS(size);
63 mp_ptr qp;
64
65 qp = MPZ_REALLOC (q, abs_size);
66
67 mpn_bdiv_dbm1 (qp, PTR(a), abs_size, GMP_NUMB_MASK / 3);
68
69 abs_size -= (qp[abs_size-1] == 0);
70 SIZ(q) = (size>0 ? abs_size : -abs_size);
71 }
72 #endif
73
74 #if GMP_NUMB_BITS % 4 == 0
75 static void
76 mpz_divexact_by5 (mpz_ptr q, mpz_srcptr a)
77 {
78 mp_size_t size = SIZ(a);
79 mp_size_t abs_size = ABS(size);
80 mp_ptr qp;
81
82 qp = MPZ_REALLOC (q, abs_size);
83
84 mpn_bdiv_dbm1 (qp, PTR(a), abs_size, GMP_NUMB_MASK / 5);
85
86 abs_size -= (qp[abs_size-1] == 0);
87 SIZ(q) = (size>0 ? abs_size : -abs_size);
88 }
89 #endif
90
91 static void
92 mpz_divexact_limb (mpz_ptr q, mpz_srcptr a, mp_limb_t d)
93 {
94 mp_size_t size = SIZ(a);
95 mp_size_t abs_size = ABS(size);
96 mp_ptr qp;
97
98 qp = MPZ_REALLOC (q, abs_size);
99
100 MPN_DIVREM_OR_DIVEXACT_1 (qp, PTR(a), abs_size, d);
101
102 abs_size -= (qp[abs_size-1] == 0);
103 SIZ(q) = (size > 0 ? abs_size : -abs_size);
104 }
105
106 void
107 mpz_divexact_gcd (mpz_ptr q, mpz_srcptr a, mpz_srcptr d)
108 {
109 ASSERT (mpz_sgn (d) > 0);
110
111 if (SIZ(a) == 0)
112 {
113 SIZ(q) = 0;
114 return;
115 }
116
117 if (SIZ(d) == 1)
118 {
119 mp_limb_t dl = PTR(d)[0];
120 int twos;
121
122 if ((dl & 1) == 0)
123 {
124 count_trailing_zeros (twos, dl);
125 dl >>= twos;
126 mpz_tdiv_q_2exp (q, a, twos);
127 a = q;
128 }
129
130 if (dl == 1)
131 {
132 if (q != a)
133 mpz_set (q, a);
134 return;
135 }
136 #if GMP_NUMB_BITS % 2 == 0
137 if (dl == 3)
138 {
139 mpz_divexact_by3 (q, a);
140 return;
141 }
142 #endif
143 #if GMP_NUMB_BITS % 4 == 0
144 if (dl == 5)
145 {
146 mpz_divexact_by5 (q, a);
147 return;
148 }
149 #endif
150
151 mpz_divexact_limb (q, a, dl);
152 return;
153 }
154
155 mpz_divexact (q, a, d);
156 }