1 /* gcd_subdiv_step.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, 2010, 2011 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 <stdlib.h> /* for NULL */
36
37 #include "gmp.h"
38 #include "gmp-impl.h"
39 #include "longlong.h"
40
41 /* Used when mpn_hgcd or mpn_hgcd2 has failed. Then either one of a or
42 b is small, or the difference is small. Perform one subtraction
43 followed by one division. The normal case is to compute the reduced
44 a and b, and return the new size.
45
46 If s == 0 (used for gcd and gcdext), returns zero if the gcd is
47 found.
48
49 If s > 0, don't reduce to size <= s, and return zero if no
50 reduction is possible (if either a, b or |a-b| is of size <= s). */
51
52 /* The hook function is called as
53
54 hook(ctx, gp, gn, qp, qn, d)
55
56 in the following cases:
57
58 + If A = B at the start, G is the gcd, Q is NULL, d = -1.
59
60 + If one input is zero at the start, G is the gcd, Q is NULL,
61 d = 0 if A = G and d = 1 if B = G.
62
63 Otherwise, if d = 0 we have just subtracted a multiple of A from B,
64 and if d = 1 we have subtracted a multiple of B from A.
65
66 + If A = B after subtraction, G is the gcd, Q is NULL.
67
68 + If we get a zero remainder after division, G is the gcd, Q is the
69 quotient.
70
71 + Otherwise, G is NULL, Q is the quotient (often 1).
72
73 */
74
75 mp_size_t
mpn_gcd_subdiv_step(mp_ptr ap,mp_ptr bp,mp_size_t n,mp_size_t s,gcd_subdiv_step_hook * hook,void * ctx,mp_ptr tp)76 mpn_gcd_subdiv_step (mp_ptr ap, mp_ptr bp, mp_size_t n, mp_size_t s,
77 gcd_subdiv_step_hook *hook, void *ctx,
78 mp_ptr tp)
79 {
80 static const mp_limb_t one = CNST_LIMB(1);
81 mp_size_t an, bn, qn;
82
83 int swapped;
84
85 ASSERT (n > 0);
86 ASSERT (ap[n-1] > 0 || bp[n-1] > 0);
87
88 an = bn = n;
89 MPN_NORMALIZE (ap, an);
90 MPN_NORMALIZE (bp, bn);
91
92 swapped = 0;
93
94 /* Arrange so that a < b, subtract b -= a, and maintain
95 normalization. */
96 if (an == bn)
97 {
98 int c;
99 MPN_CMP (c, ap, bp, an);
100 if (UNLIKELY (c == 0))
101 {
102 /* For gcdext, return the smallest of the two cofactors, so
103 pass d = -1. */
104 if (s == 0)
105 hook (ctx, ap, an, NULL, 0, -1);
106 return 0;
107 }
108 else if (c > 0)
109 {
110 MP_PTR_SWAP (ap, bp);
111 swapped ^= 1;
112 }
113 }
114 else
115 {
116 if (an > bn)
117 {
118 MPN_PTR_SWAP (ap, an, bp, bn);
119 swapped ^= 1;
120 }
121 }
122 if (an <= s)
123 {
124 if (s == 0)
125 hook (ctx, bp, bn, NULL, 0, swapped ^ 1);
126 return 0;
127 }
128
129 ASSERT_NOCARRY (mpn_sub (bp, bp, bn, ap, an));
130 MPN_NORMALIZE (bp, bn);
131 ASSERT (bn > 0);
132
133 if (bn <= s)
134 {
135 /* Undo subtraction. */
136 mp_limb_t cy = mpn_add (bp, ap, an, bp, bn);
137 if (cy > 0)
138 bp[an] = cy;
139 return 0;
140 }
141
142 /* Arrange so that a < b */
143 if (an == bn)
144 {
145 int c;
146 MPN_CMP (c, ap, bp, an);
147 if (UNLIKELY (c == 0))
148 {
149 if (s > 0)
150 /* Just record subtraction and return */
151 hook (ctx, NULL, 0, &one, 1, swapped);
152 else
153 /* Found gcd. */
154 hook (ctx, bp, bn, NULL, 0, swapped);
155 return 0;
156 }
157
158 hook (ctx, NULL, 0, &one, 1, swapped);
159
160 if (c > 0)
161 {
162 MP_PTR_SWAP (ap, bp);
163 swapped ^= 1;
164 }
165 }
166 else
167 {
168 hook (ctx, NULL, 0, &one, 1, swapped);
169
170 if (an > bn)
171 {
172 MPN_PTR_SWAP (ap, an, bp, bn);
173 swapped ^= 1;
174 }
175 }
176
177 mpn_tdiv_qr (tp, bp, 0, bp, bn, ap, an);
178 qn = bn - an + 1;
179 bn = an;
180 MPN_NORMALIZE (bp, bn);
181
182 if (UNLIKELY (bn <= s))
183 {
184 if (s == 0)
185 {
186 hook (ctx, ap, an, tp, qn, swapped);
187 return 0;
188 }
189
190 /* Quotient is one too large, so decrement it and add back A. */
191 if (bn > 0)
192 {
193 mp_limb_t cy = mpn_add (bp, ap, an, bp, bn);
194 if (cy)
195 bp[an++] = cy;
196 }
197 else
198 MPN_COPY (bp, ap, an);
199
200 MPN_DECR_U (tp, qn, 1);
201 }
202
203 hook (ctx, NULL, 0, tp, qn, swapped);
204 return an;
205 }
206