1 /* hgcd.c.
2
3 THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES. IT IS ONLY
4 SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
5 GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
6
7 Copyright 2003-2005, 2008, 2011, 2012 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 either:
13
14 * the GNU Lesser General Public License as published by the Free
15 Software Foundation; either version 3 of the License, or (at your
16 option) any later version.
17
18 or
19
20 * the GNU General Public License as published by the Free Software
21 Foundation; either version 2 of the License, or (at your option) any
22 later version.
23
24 or both in parallel, as here.
25
26 The GNU MP Library is distributed in the hope that it will be useful, but
27 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
28 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
29 for more details.
30
31 You should have received copies of the GNU General Public License and the
32 GNU Lesser General Public License along with the GNU MP Library. If not,
33 see https://www.gnu.org/licenses/. */
34
35 #include "gmp-impl.h"
36 #include "longlong.h"
37
38
39 /* Size analysis for hgcd:
40
41 For the recursive calls, we have n1 <= ceil(n / 2). Then the
42 storage need is determined by the storage for the recursive call
43 computing M1, and hgcd_matrix_adjust and hgcd_matrix_mul calls that use M1
44 (after this, the storage needed for M1 can be recycled).
45
46 Let S(r) denote the required storage. For M1 we need 4 * (ceil(n1/2) + 1)
47 = 4 * (ceil(n/4) + 1), for the hgcd_matrix_adjust call, we need n + 2,
48 and for the hgcd_matrix_mul, we may need 3 ceil(n/2) + 8. In total,
49 4 * ceil(n/4) + 3 ceil(n/2) + 12 <= 10 ceil(n/4) + 12.
50
51 For the recursive call, we need S(n1) = S(ceil(n/2)).
52
53 S(n) <= 10*ceil(n/4) + 12 + S(ceil(n/2))
54 <= 10*(ceil(n/4) + ... + ceil(n/2^(1+k))) + 12k + S(ceil(n/2^k))
55 <= 10*(2 ceil(n/4) + k) + 12k + S(ceil(n/2^k))
56 <= 20 ceil(n/4) + 22k + S(ceil(n/2^k))
57 */
58
59 mp_size_t
60 mpn_hgcd_itch (mp_size_t n)
61 {
62 unsigned k;
63 int count;
64 mp_size_t nscaled;
65
66 if (BELOW_THRESHOLD (n, HGCD_THRESHOLD))
67 return n;
68
69 /* Get the recursion depth. */
70 nscaled = (n - 1) / (HGCD_THRESHOLD - 1);
71 count_leading_zeros (count, nscaled);
72 k = GMP_LIMB_BITS - count;
73
74 return 20 * ((n+3) / 4) + 22 * k + HGCD_THRESHOLD;
75 }
76
77 /* Reduces a,b until |a-b| fits in n/2 + 1 limbs. Constructs matrix M
78 with elements of size at most (n+1)/2 - 1. Returns new size of a,
79 b, or zero if no reduction is possible. */
80
81 mp_size_t
82 mpn_hgcd (mp_ptr ap, mp_ptr bp, mp_size_t n,
83 struct hgcd_matrix *M, mp_ptr tp)
84 {
85 mp_size_t s = n/2 + 1;
86
87 mp_size_t nn;
88 int success = 0;
89
90 if (n <= s)
91 /* Happens when n <= 2, a fairly uninteresting case but exercised
92 by the random inputs of the testsuite. */
93 return 0;
94
95 ASSERT ((ap[n-1] | bp[n-1]) > 0);
96
97 ASSERT ((n+1)/2 - 1 < M->alloc);
98
99 if (ABOVE_THRESHOLD (n, HGCD_THRESHOLD))
100 {
101 mp_size_t n2 = (3*n)/4 + 1;
102 mp_size_t p = n/2;
103
104 nn = mpn_hgcd_reduce (M, ap, bp, n, p, tp);
105 if (nn)
106 {
107 n = nn;
108 success = 1;
109 }
110
111 /* NOTE: It appears this loop never runs more than once (at
112 least when not recursing to hgcd_appr). */
113 while (n > n2)
114 {
115 /* Needs n + 1 storage */
116 nn = mpn_hgcd_step (n, ap, bp, s, M, tp);
117 if (!nn)
118 return success ? n : 0;
119
120 n = nn;
121 success = 1;
122 }
123
124 if (n > s + 2)
125 {
126 struct hgcd_matrix M1;
127 mp_size_t scratch;
128
129 p = 2*s - n + 1;
130 scratch = MPN_HGCD_MATRIX_INIT_ITCH (n-p);
131
132 mpn_hgcd_matrix_init(&M1, n - p, tp);
133
134 /* FIXME: Should use hgcd_reduce, but that may require more
135 scratch space, which requires review. */
136
137 nn = mpn_hgcd (ap + p, bp + p, n - p, &M1, tp + scratch);
138 if (nn > 0)
139 {
140 /* We always have max(M) > 2^{-(GMP_NUMB_BITS + 1)} max(M1) */
141 ASSERT (M->n + 2 >= M1.n);
142
143 /* Furthermore, assume M ends with a quotient (1, q; 0, 1),
144 then either q or q + 1 is a correct quotient, and M1 will
145 start with either (1, 0; 1, 1) or (2, 1; 1, 1). This
146 rules out the case that the size of M * M1 is much
147 smaller than the expected M->n + M1->n. */
148
149 ASSERT (M->n + M1.n < M->alloc);
150
151 /* Needs 2 (p + M->n) <= 2 (2*s - n2 + 1 + n2 - s - 1)
152 = 2*s <= 2*(floor(n/2) + 1) <= n + 2. */
153 n = mpn_hgcd_matrix_adjust (&M1, p + nn, ap, bp, p, tp + scratch);
154
155 /* We need a bound for of M->n + M1.n. Let n be the original
156 input size. Then
157
158 ceil(n/2) - 1 >= size of product >= M.n + M1.n - 2
159
160 and it follows that
161
162 M.n + M1.n <= ceil(n/2) + 1
163
164 Then 3*(M.n + M1.n) + 5 <= 3 * ceil(n/2) + 8 is the
165 amount of needed scratch space. */
166 mpn_hgcd_matrix_mul (M, &M1, tp + scratch);
167 success = 1;
168 }
169 }
170 }
171
172 for (;;)
173 {
174 /* Needs s+3 < n */
175 nn = mpn_hgcd_step (n, ap, bp, s, M, tp);
176 if (!nn)
177 return success ? n : 0;
178
179 n = nn;
180 success = 1;
181 }
182 }