1 /* gcdext_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, 2004, 2005, 2008, 2009 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 the GNU Lesser General Public License as published by
13 the Free Software Foundation; either version 3 of the License, or (at your
14 option) any later version.
15
16 The GNU MP Library is distributed in the hope that it will be useful, but
17 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19 License for more details.
20
21 You should have received a copy of the GNU Lesser General Public License
22 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
23
24 #include "gmp.h"
25 #include "gmp-impl.h"
26 #include "longlong.h"
27
28 /* Used when mpn_hgcd or mpn_hgcd2 has failed. Then either one of a or
29 b is small, or the difference is small. Perform one subtraction
30 followed by one division. If the gcd is found, stores it in gp and
31 *gn, and returns zero. Otherwise, compute the reduced a and b,
32 return the new size, and cofactors. */
33
34 /* Temporary storage: Needs n limbs for the quotient, at qp. tp must
35 point to an area large enough for the resulting cofactor, plus one
36 limb extra. All in all, 2N + 1 if N is a bound for both inputs and
37 outputs. */
38 mp_size_t
mpn_gcdext_subdiv_step(mp_ptr gp,mp_size_t * gn,mp_ptr up,mp_size_t * usizep,mp_ptr ap,mp_ptr bp,mp_size_t n,mp_ptr u0,mp_ptr u1,mp_size_t * unp,mp_ptr qp,mp_ptr tp)39 mpn_gcdext_subdiv_step (mp_ptr gp, mp_size_t *gn, mp_ptr up, mp_size_t *usizep,
40 mp_ptr ap, mp_ptr bp, mp_size_t n,
41 mp_ptr u0, mp_ptr u1, mp_size_t *unp,
42 mp_ptr qp, mp_ptr tp)
43 {
44 mp_size_t an, bn, un;
45 mp_size_t qn;
46 mp_size_t u0n;
47
48 int swapped;
49
50 an = bn = n;
51
52 ASSERT (an > 0);
53 ASSERT (ap[an-1] > 0 || bp[an-1] > 0);
54
55 MPN_NORMALIZE (ap, an);
56 MPN_NORMALIZE (bp, bn);
57
58 un = *unp;
59
60 swapped = 0;
61
62 if (UNLIKELY (an == 0))
63 {
64 return_b:
65 MPN_COPY (gp, bp, bn);
66 *gn = bn;
67
68 MPN_NORMALIZE (u0, un);
69 MPN_COPY (up, u0, un);
70
71 *usizep = swapped ? un : -un;
72
73 return 0;
74 }
75 else if (UNLIKELY (bn == 0))
76 {
77 MPN_COPY (gp, ap, an);
78 *gn = an;
79
80 MPN_NORMALIZE (u1, un);
81 MPN_COPY (up, u1, un);
82
83 *usizep = swapped ? -un : un;
84
85 return 0;
86 }
87
88 /* Arrange so that a > b, subtract an -= bn, and maintain
89 normalization. */
90 if (an < bn)
91 {
92 MPN_PTR_SWAP (ap, an, bp, bn);
93 MP_PTR_SWAP (u0, u1);
94 swapped ^= 1;
95 }
96 else if (an == bn)
97 {
98 int c;
99 MPN_CMP (c, ap, bp, an);
100 if (UNLIKELY (c == 0))
101 {
102 MPN_COPY (gp, ap, an);
103 *gn = an;
104
105 /* Must return the smallest cofactor, +u1 or -u0 */
106 MPN_CMP (c, u0, u1, un);
107 ASSERT (c != 0 || (un == 1 && u0[0] == 1 && u1[0] == 1));
108
109 if (c < 0)
110 {
111 MPN_NORMALIZE (u0, un);
112 MPN_COPY (up, u0, un);
113 swapped ^= 1;
114 }
115 else
116 {
117 MPN_NORMALIZE_NOT_ZERO (u1, un);
118 MPN_COPY (up, u1, un);
119 }
120
121 *usizep = swapped ? -un : un;
122 return 0;
123 }
124 else if (c < 0)
125 {
126 MP_PTR_SWAP (ap, bp);
127 MP_PTR_SWAP (u0, u1);
128 swapped ^= 1;
129 }
130 }
131 /* Reduce a -= b, u1 += u0 */
132 ASSERT_NOCARRY (mpn_sub (ap, ap, an, bp, bn));
133 MPN_NORMALIZE (ap, an);
134 ASSERT (an > 0);
135
136 u1[un] = mpn_add_n (u1, u1, u0, un);
137 un += (u1[un] > 0);
138
139 /* Arrange so that a > b, and divide a = q b + r */
140 if (an < bn)
141 {
142 MPN_PTR_SWAP (ap, an, bp, bn);
143 MP_PTR_SWAP (u0, u1);
144 swapped ^= 1;
145 }
146 else if (an == bn)
147 {
148 int c;
149 MPN_CMP (c, ap, bp, an);
150 if (UNLIKELY (c == 0))
151 goto return_b;
152 else if (c < 0)
153 {
154 MP_PTR_SWAP (ap, bp);
155 MP_PTR_SWAP (u0, u1);
156 swapped ^= 1;
157 }
158 }
159
160 /* Reduce a -= q b, u1 += q u0 */
161 qn = an - bn + 1;
162 mpn_tdiv_qr (qp, ap, 0, ap, an, bp, bn);
163
164 if (mpn_zero_p (ap, bn))
165 goto return_b;
166
167 n = bn;
168
169 /* Update u1 += q u0 */
170 u0n = un;
171 MPN_NORMALIZE (u0, u0n);
172
173 if (u0n > 0)
174 {
175 qn -= (qp[qn - 1] == 0);
176
177 if (qn > u0n)
178 mpn_mul (tp, qp, qn, u0, u0n);
179 else
180 mpn_mul (tp, u0, u0n, qp, qn);
181
182 if (qn + u0n > un)
183 {
184 mp_size_t u1n = un;
185 un = qn + u0n;
186 un -= (tp[un-1] == 0);
187 u1[un] = mpn_add (u1, tp, un, u1, u1n);
188 }
189 else
190 {
191 u1[un] = mpn_add (u1, u1, un, tp, qn + u0n);
192 }
193
194 un += (u1[un] > 0);
195 }
196
197 *unp = un;
198 return n;
199 }
200