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