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; | |
3d6115d1 | 54 | mp *q = MP_NEW, *t, *spare = MP_NEW; |
ceb3f0c0 | 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)) { | |
3d6115d1 | 72 | t = a; a = b; b = t; |
ceb3f0c0 | 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 | gf_div(&q, &u, u, v); |
116 | if (f & f_ext) { | |
3d6115d1 | 117 | t = gf_mul(spare, X, q); |
cf76bcbb | 118 | t = gf_add(t, t, x); |
3d6115d1 MW |
119 | spare = x; x = X; X = t; |
120 | t = gf_mul(spare, Y, q); | |
cf76bcbb | 121 | t = gf_add(t, t, y); |
3d6115d1 | 122 | spare = y; y = Y; Y = t; |
ceb3f0c0 | 123 | } |
124 | t = u; u = v; v = t; | |
125 | } | |
126 | ||
3d6115d1 | 127 | MP_DROP(q); if (spare) MP_DROP(spare); |
ceb3f0c0 | 128 | if (!gcd) |
129 | MP_DROP(u); | |
130 | else { | |
131 | if (*gcd) MP_DROP(*gcd); | |
132 | u->f &= ~MP_NEG; | |
133 | *gcd = u; | |
134 | } | |
135 | ||
136 | /* --- Perform a little normalization --- */ | |
137 | ||
138 | if (f & f_ext) { | |
139 | ||
140 | /* --- If @a@ and @b@ got swapped, swap the coefficients back --- */ | |
141 | ||
142 | if (f & f_swap) { | |
3d6115d1 | 143 | t = x; x = y; y = t; |
ceb3f0c0 | 144 | t = a; a = b; b = t; |
145 | } | |
146 | ||
147 | /* --- Store the results --- */ | |
148 | ||
149 | if (!xx) | |
150 | MP_DROP(x); | |
151 | else { | |
152 | if (*xx) MP_DROP(*xx); | |
153 | *xx = x; | |
154 | } | |
155 | ||
156 | if (!yy) | |
157 | MP_DROP(y); | |
158 | else { | |
159 | if (*yy) MP_DROP(*yy); | |
160 | *yy = y; | |
161 | } | |
162 | } | |
163 | ||
164 | MP_DROP(v); | |
165 | MP_DROP(X); MP_DROP(Y); | |
ceb3f0c0 | 166 | } |
167 | ||
b817bfc6 | 168 | /* -- @gf_modinv@ --- * |
169 | * | |
170 | * Arguments: @mp *d@ = destination | |
171 | * @mp *x@ = argument | |
172 | * @mp *p@ = modulus | |
173 | * | |
174 | * Returns: The inverse %$x^{-1} \bmod p$%. | |
175 | * | |
176 | * Use: Computes a modular inverse, the catch being that the | |
177 | * arguments and results are binary polynomials. An assertion | |
178 | * fails if %$p$% has no inverse. | |
179 | */ | |
180 | ||
181 | mp *gf_modinv(mp *d, mp *x, mp *p) | |
182 | { | |
183 | mp *g = MP_NEW; | |
184 | gf_gcd(&g, 0, &d, p, x); | |
185 | assert(MP_EQ(g, MP_ONE)); | |
186 | mp_drop(g); | |
187 | return (d); | |
188 | } | |
189 | ||
ceb3f0c0 | 190 | /*----- Test rig ----------------------------------------------------------*/ |
191 | ||
192 | #ifdef TEST_RIG | |
193 | ||
194 | static int gcd(dstr *v) | |
195 | { | |
196 | int ok = 1; | |
197 | mp *a = *(mp **)v[0].buf; | |
198 | mp *b = *(mp **)v[1].buf; | |
199 | mp *g = *(mp **)v[2].buf; | |
200 | mp *x = *(mp **)v[3].buf; | |
201 | mp *y = *(mp **)v[4].buf; | |
202 | ||
203 | mp *gg = MP_NEW, *xx = MP_NEW, *yy = MP_NEW; | |
204 | gf_gcd(&gg, &xx, &yy, a, b); | |
205 | if (!MP_EQ(x, xx)) { | |
b817bfc6 | 206 | fputs("\n*** gf_gcd(x) failed", stderr); |
45c0fd36 MW |
207 | fputs("\na = ", stderr); mp_writefile(a, stderr, 16); |
208 | fputs("\nb = ", stderr); mp_writefile(b, stderr, 16); | |
ceb3f0c0 | 209 | fputs("\nexpect = ", stderr); mp_writefile(x, stderr, 16); |
210 | fputs("\nresult = ", stderr); mp_writefile(xx, stderr, 16); | |
211 | fputc('\n', stderr); | |
212 | ok = 0; | |
213 | } | |
214 | if (!MP_EQ(y, yy)) { | |
b817bfc6 | 215 | fputs("\n*** gf_gcd(y) failed", stderr); |
45c0fd36 MW |
216 | fputs("\na = ", stderr); mp_writefile(a, stderr, 16); |
217 | fputs("\nb = ", stderr); mp_writefile(b, stderr, 16); | |
ceb3f0c0 | 218 | fputs("\nexpect = ", stderr); mp_writefile(y, stderr, 16); |
219 | fputs("\nresult = ", stderr); mp_writefile(yy, stderr, 16); | |
220 | fputc('\n', stderr); | |
221 | ok = 0; | |
222 | } | |
223 | ||
224 | if (!ok) { | |
225 | mp *ax = gf_mul(MP_NEW, a, xx); | |
226 | mp *by = gf_mul(MP_NEW, b, yy); | |
227 | ax = gf_add(ax, ax, by); | |
228 | if (MP_EQ(ax, gg)) | |
229 | fputs("\n*** (Alternative result found.)\n", stderr); | |
230 | MP_DROP(ax); | |
231 | MP_DROP(by); | |
232 | } | |
233 | ||
234 | if (!MP_EQ(g, gg)) { | |
b817bfc6 | 235 | fputs("\n*** gf_gcd(gcd) failed", stderr); |
45c0fd36 MW |
236 | fputs("\na = ", stderr); mp_writefile(a, stderr, 16); |
237 | fputs("\nb = ", stderr); mp_writefile(b, stderr, 16); | |
ceb3f0c0 | 238 | fputs("\nexpect = ", stderr); mp_writefile(g, stderr, 16); |
239 | fputs("\nresult = ", stderr); mp_writefile(gg, stderr, 16); | |
240 | fputc('\n', stderr); | |
241 | ok = 0; | |
242 | } | |
243 | MP_DROP(a); MP_DROP(b); MP_DROP(g); MP_DROP(x); MP_DROP(y); | |
244 | MP_DROP(gg); MP_DROP(xx); MP_DROP(yy); | |
245 | assert(mparena_count(MPARENA_GLOBAL) == 0); | |
246 | return (ok); | |
247 | } | |
248 | ||
249 | static test_chunk tests[] = { | |
250 | { "gcd", gcd, { &type_mp, &type_mp, &type_mp, &type_mp, &type_mp, 0 } }, | |
251 | { 0, 0, { 0 } } | |
252 | }; | |
253 | ||
254 | int main(int argc, char *argv[]) | |
255 | { | |
256 | sub_init(); | |
257 | test_run(argc, argv, tests, SRCDIR "/tests/gf"); | |
258 | return (0); | |
259 | } | |
260 | ||
261 | #endif | |
262 | ||
263 | /*----- That's all, folks -------------------------------------------------*/ |