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