Commit | Line | Data |
---|---|---|
ceb3f0c0 | 1 | /* -*-c-*- |
2 | * | |
a69a3efd | 3 | * $Id$ |
ceb3f0c0 | 4 | * |
5 | * Euclidian algorithm on binary polynomials | |
6 | * | |
7 | * (c) 2004 Straylight/Edgeware | |
8 | */ | |
9 | ||
45c0fd36 | 10 | /*----- Licensing notice --------------------------------------------------* |
ceb3f0c0 | 11 | * |
12 | * This file is part of Catacomb. | |
13 | * | |
14 | * Catacomb is free software; you can redistribute it and/or modify | |
15 | * it under the terms of the GNU Library General Public License as | |
16 | * published by the Free Software Foundation; either version 2 of the | |
17 | * License, or (at your option) any later version. | |
45c0fd36 | 18 | * |
ceb3f0c0 | 19 | * Catacomb is distributed in the hope that it will be useful, |
20 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | |
21 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
22 | * GNU Library General Public License for more details. | |
45c0fd36 | 23 | * |
ceb3f0c0 | 24 | * You should have received a copy of the GNU Library General Public |
25 | * License along with Catacomb; if not, write to the Free | |
26 | * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, | |
27 | * MA 02111-1307, USA. | |
28 | */ | |
29 | ||
ceb3f0c0 | 30 | /*----- Header files ------------------------------------------------------*/ |
31 | ||
32 | #include "gf.h" | |
33 | ||
34 | /*----- Main code ---------------------------------------------------------*/ | |
35 | ||
36 | /* --- @gf_gcd@ --- * | |
37 | * | |
38 | * Arguments: @mp **gcd, **xx, **yy@ = where to write the results | |
39 | * @mp *a, *b@ = sources (must be nonzero) | |
40 | * | |
41 | * | |
42 | * Returns: --- | |
43 | * | |
44 | * Use: Calculates @gcd(a, b)@, and two numbers @x@ and @y@ such that | |
45 | * @ax + by = gcd(a, b)@. This is useful for computing modular | |
46 | * inverses. | |
47 | */ | |
48 | ||
49 | void gf_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b) | |
50 | { | |
51 | mp *x = MP_ONE, *X = MP_ZERO; | |
52 | mp *y = MP_ZERO, *Y = MP_ONE; | |
53 | mp *u, *v; | |
54 | mp *q = MP_NEW; | |
55 | unsigned f = 0; | |
56 | ||
57 | #define f_swap 1u | |
58 | #define f_ext 2u | |
59 | ||
60 | /* --- Sort out some initial flags --- */ | |
61 | ||
62 | if (xx || yy) | |
63 | f |= f_ext; | |
64 | ||
65 | /* --- Ensure that @a@ is larger than @b@ --- * | |
66 | * | |
67 | * If they're the same length we don't care which order they're in, so this | |
68 | * unsigned comparison is fine. | |
69 | */ | |
70 | ||
71 | if (MPX_UCMP(a->v, a->vl, <, b->v, b->vl)) { | |
72 | { mp *t = a; a = b; b = t; } | |
73 | f |= f_swap; | |
74 | } | |
75 | ||
76 | /* --- Check for zeroness --- */ | |
77 | ||
78 | if (MP_EQ(b, MP_ZERO)) { | |
79 | ||
80 | /* --- Store %$|a|$% as the GCD --- */ | |
81 | ||
82 | if (gcd) { | |
83 | if (*gcd) MP_DROP(*gcd); | |
84 | a = MP_COPY(a); | |
85 | *gcd = a; | |
86 | } | |
87 | ||
88 | /* --- Store %$1$% and %$0$% in the appropriate bins --- */ | |
89 | ||
90 | if (f & f_ext) { | |
91 | if (f & f_swap) { | |
92 | mp **t = xx; xx = yy; yy = t; | |
93 | } | |
94 | if (xx) { | |
95 | if (*xx) MP_DROP(*xx); | |
96 | if (MP_EQ(a, MP_ZERO)) | |
97 | *xx = MP_ZERO; | |
98 | else | |
99 | *xx = MP_ONE; | |
100 | } | |
101 | if (yy) { | |
102 | if (*yy) MP_DROP(*yy); | |
103 | *yy = MP_ZERO; | |
104 | } | |
105 | } | |
106 | return; | |
107 | } | |
108 | ||
cf76bcbb | 109 | /* --- Main extended Euclidean algorithm --- */ |
ceb3f0c0 | 110 | |
111 | u = MP_COPY(a); | |
112 | v = MP_COPY(b); | |
113 | ||
a69a3efd | 114 | while (!MP_ZEROP(v)) { |
ceb3f0c0 | 115 | mp *t; |
116 | gf_div(&q, &u, u, v); | |
117 | if (f & f_ext) { | |
118 | t = gf_mul(MP_NEW, X, q); | |
cf76bcbb | 119 | t = gf_add(t, t, x); |
ceb3f0c0 | 120 | MP_DROP(x); x = X; X = t; |
121 | t = gf_mul(MP_NEW, Y, q); | |
cf76bcbb | 122 | t = gf_add(t, t, y); |
ceb3f0c0 | 123 | MP_DROP(y); y = Y; Y = t; |
124 | } | |
125 | t = u; u = v; v = t; | |
126 | } | |
127 | ||
128 | MP_DROP(q); | |
129 | if (!gcd) | |
130 | MP_DROP(u); | |
131 | else { | |
132 | if (*gcd) MP_DROP(*gcd); | |
133 | u->f &= ~MP_NEG; | |
134 | *gcd = u; | |
135 | } | |
136 | ||
137 | /* --- Perform a little normalization --- */ | |
138 | ||
139 | if (f & f_ext) { | |
140 | ||
141 | /* --- If @a@ and @b@ got swapped, swap the coefficients back --- */ | |
142 | ||
143 | if (f & f_swap) { | |
144 | mp *t = x; x = y; y = t; | |
145 | t = a; a = b; b = t; | |
146 | } | |
147 | ||
148 | /* --- Store the results --- */ | |
149 | ||
150 | if (!xx) | |
151 | MP_DROP(x); | |
152 | else { | |
153 | if (*xx) MP_DROP(*xx); | |
154 | *xx = x; | |
155 | } | |
156 | ||
157 | if (!yy) | |
158 | MP_DROP(y); | |
159 | else { | |
160 | if (*yy) MP_DROP(*yy); | |
161 | *yy = y; | |
162 | } | |
163 | } | |
164 | ||
165 | MP_DROP(v); | |
166 | MP_DROP(X); MP_DROP(Y); | |
ceb3f0c0 | 167 | } |
168 | ||
b817bfc6 | 169 | /* -- @gf_modinv@ --- * |
170 | * | |
171 | * Arguments: @mp *d@ = destination | |
172 | * @mp *x@ = argument | |
173 | * @mp *p@ = modulus | |
174 | * | |
175 | * Returns: The inverse %$x^{-1} \bmod p$%. | |
176 | * | |
177 | * Use: Computes a modular inverse, the catch being that the | |
178 | * arguments and results are binary polynomials. An assertion | |
179 | * fails if %$p$% has no inverse. | |
180 | */ | |
181 | ||
182 | mp *gf_modinv(mp *d, mp *x, mp *p) | |
183 | { | |
184 | mp *g = MP_NEW; | |
185 | gf_gcd(&g, 0, &d, p, x); | |
186 | assert(MP_EQ(g, MP_ONE)); | |
187 | mp_drop(g); | |
188 | return (d); | |
189 | } | |
190 | ||
ceb3f0c0 | 191 | /*----- Test rig ----------------------------------------------------------*/ |
192 | ||
193 | #ifdef TEST_RIG | |
194 | ||
195 | static int gcd(dstr *v) | |
196 | { | |
197 | int ok = 1; | |
198 | mp *a = *(mp **)v[0].buf; | |
199 | mp *b = *(mp **)v[1].buf; | |
200 | mp *g = *(mp **)v[2].buf; | |
201 | mp *x = *(mp **)v[3].buf; | |
202 | mp *y = *(mp **)v[4].buf; | |
203 | ||
204 | mp *gg = MP_NEW, *xx = MP_NEW, *yy = MP_NEW; | |
205 | gf_gcd(&gg, &xx, &yy, a, b); | |
206 | if (!MP_EQ(x, xx)) { | |
b817bfc6 | 207 | fputs("\n*** gf_gcd(x) failed", stderr); |
45c0fd36 MW |
208 | fputs("\na = ", stderr); mp_writefile(a, stderr, 16); |
209 | fputs("\nb = ", stderr); mp_writefile(b, stderr, 16); | |
ceb3f0c0 | 210 | fputs("\nexpect = ", stderr); mp_writefile(x, stderr, 16); |
211 | fputs("\nresult = ", stderr); mp_writefile(xx, stderr, 16); | |
212 | fputc('\n', stderr); | |
213 | ok = 0; | |
214 | } | |
215 | if (!MP_EQ(y, yy)) { | |
b817bfc6 | 216 | fputs("\n*** gf_gcd(y) failed", stderr); |
45c0fd36 MW |
217 | fputs("\na = ", stderr); mp_writefile(a, stderr, 16); |
218 | fputs("\nb = ", stderr); mp_writefile(b, stderr, 16); | |
ceb3f0c0 | 219 | fputs("\nexpect = ", stderr); mp_writefile(y, stderr, 16); |
220 | fputs("\nresult = ", stderr); mp_writefile(yy, stderr, 16); | |
221 | fputc('\n', stderr); | |
222 | ok = 0; | |
223 | } | |
224 | ||
225 | if (!ok) { | |
226 | mp *ax = gf_mul(MP_NEW, a, xx); | |
227 | mp *by = gf_mul(MP_NEW, b, yy); | |
228 | ax = gf_add(ax, ax, by); | |
229 | if (MP_EQ(ax, gg)) | |
230 | fputs("\n*** (Alternative result found.)\n", stderr); | |
231 | MP_DROP(ax); | |
232 | MP_DROP(by); | |
233 | } | |
234 | ||
235 | if (!MP_EQ(g, gg)) { | |
b817bfc6 | 236 | fputs("\n*** gf_gcd(gcd) failed", stderr); |
45c0fd36 MW |
237 | fputs("\na = ", stderr); mp_writefile(a, stderr, 16); |
238 | fputs("\nb = ", stderr); mp_writefile(b, stderr, 16); | |
ceb3f0c0 | 239 | fputs("\nexpect = ", stderr); mp_writefile(g, stderr, 16); |
240 | fputs("\nresult = ", stderr); mp_writefile(gg, stderr, 16); | |
241 | fputc('\n', stderr); | |
242 | ok = 0; | |
243 | } | |
244 | MP_DROP(a); MP_DROP(b); MP_DROP(g); MP_DROP(x); MP_DROP(y); | |
245 | MP_DROP(gg); MP_DROP(xx); MP_DROP(yy); | |
246 | assert(mparena_count(MPARENA_GLOBAL) == 0); | |
247 | return (ok); | |
248 | } | |
249 | ||
250 | static test_chunk tests[] = { | |
251 | { "gcd", gcd, { &type_mp, &type_mp, &type_mp, &type_mp, &type_mp, 0 } }, | |
252 | { 0, 0, { 0 } } | |
253 | }; | |
254 | ||
255 | int main(int argc, char *argv[]) | |
256 | { | |
257 | sub_init(); | |
258 | test_run(argc, argv, tests, SRCDIR "/tests/gf"); | |
259 | return (0); | |
260 | } | |
261 | ||
262 | #endif | |
263 | ||
264 | /*----- That's all, folks -------------------------------------------------*/ |