Commit | Line | Data |
---|---|---|
9f11b970 | 1 | /* -*-c-*- |
9f11b970 | 2 | * |
3 | * Compute square roots modulo a prime | |
4 | * | |
5 | * (c) 2000 Straylight/Edgeware | |
6 | */ | |
7 | ||
45c0fd36 | 8 | /*----- Licensing notice --------------------------------------------------* |
9f11b970 | 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 | * |
9f11b970 | 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 | * |
9f11b970 | 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 | ||
9f11b970 | 28 | /*----- Header files ------------------------------------------------------*/ |
29 | ||
30 | #include "fibrand.h" | |
31 | #include "grand.h" | |
32 | #include "mp.h" | |
33 | #include "mpmont.h" | |
34 | #include "mprand.h" | |
35 | ||
36 | /*----- Main code ---------------------------------------------------------*/ | |
37 | ||
38 | /* --- @mp_modsqrt@ --- * | |
39 | * | |
40 | * Arguments: @mp *d@ = destination integer | |
41 | * @mp *a@ = source integer | |
42 | * @mp *p@ = modulus (must be prime) | |
43 | * | |
44 | * Returns: If %$a$% is a quadratic residue, a square root of %$a$%; else | |
45 | * a null pointer. | |
46 | * | |
47 | * Use: Returns an integer %$x$% such that %$x^2 \equiv a \pmod{p}$%, | |
48 | * if one exists; else a null pointer. This function will not | |
49 | * work if %$p$% is composite: you must factor the modulus, take | |
50 | * a square root mod each factor, and recombine the results | |
51 | * using the Chinese Remainder Theorem. | |
222c8a43 MW |
52 | * |
53 | * We guarantee that the square root returned is the smallest | |
54 | * one (i.e., the `positive' square root). | |
9f11b970 | 55 | */ |
56 | ||
57 | mp *mp_modsqrt(mp *d, mp *a, mp *p) | |
58 | { | |
59 | mpmont mm; | |
5b295848 MW |
60 | size_t i, s; |
61 | mp *b, *c; | |
9f11b970 | 62 | mp *ainv; |
5b295848 MW |
63 | mp *r, *A, *aa; |
64 | mp *t; | |
65 | grand *gr; | |
d032072f | 66 | int j; |
9f11b970 | 67 | |
68 | /* --- Cope if %$a \not\in Q_p$% --- */ | |
69 | ||
d032072f MW |
70 | j = mp_jacobi(a, p); |
71 | if (j == -1) { | |
f1140c41 | 72 | mp_drop(d); |
9f11b970 | 73 | return (0); |
d032072f MW |
74 | } else if (j == 0) { |
75 | if (d != a) mp_drop(d); | |
76 | d = MP_COPY(a); | |
77 | return (d); | |
9f11b970 | 78 | } |
79 | ||
80 | /* --- Choose some quadratic non-residue --- */ | |
81 | ||
5b295848 MW |
82 | gr = fibrand_create(0); |
83 | b = MP_NEW; | |
84 | do b = mprand_range(b, p, gr, 0); while (mp_jacobi(b, p) != -1); | |
85 | gr->ops->destroy(gr); | |
45c0fd36 | 86 | |
5b295848 | 87 | /* --- Some initial setup --- */ |
9f11b970 | 88 | |
89 | mpmont_create(&mm, p); | |
5b295848 MW |
90 | ainv = mp_modinv(MP_NEW, a, p); /* %$a^{-1} \bmod p$% */ |
91 | ainv = mpmont_mul(&mm, ainv, ainv, mm.r2); | |
92 | t = mp_sub(MP_NEW, p, MP_ONE); | |
93 | t = mp_odd(t, t, &s); /* %$2^s t = p - 1$% */ | |
b0b682aa | 94 | b = mpmont_mul(&mm, b, b, mm.r2); |
5b295848 | 95 | c = mpmont_expr(&mm, b, b, t); /* %$b^t \bmod p$% */ |
9f11b970 | 96 | t = mp_add(t, t, MP_ONE); |
5b295848 MW |
97 | t = mp_lsr(t, t, 1); /* %$(t + 1)/2$% */ |
98 | a = mpmont_mul(&mm, MP_NEW, a, mm.r2); | |
99 | r = mpmont_expr(&mm, a, a, t); /* %$a^{(t+1)/2} \bmod p$% */ | |
9f11b970 | 100 | |
5b295848 MW |
101 | /* --- Now for the main loop --- * |
102 | * | |
103 | * Let %$g = c^{-2}$%; we know that %$g$% is a generator of the order- | |
104 | * %$2^{s-1}$% subgroup mod %$p$%. We also know that %$A = a^t = r^2/a$% | |
105 | * is an element of this group. If we can determine %$m$% such that | |
106 | * %$g^m = A$% then %$a^{(t+1)/2}/g^{m/2} = r c^m$% is the square root we | |
107 | * seek. | |
108 | * | |
109 | * Write %$m = m_0 + 2 m'$%. Then %$A^{2^{s-1}} = g^{m_0 2^{s-1}}$%, which | |
110 | * is %$1$% if %$m_0 = 0$% or %$-1$% if %$m_0 = 1$% (modulo %$p$%). Then | |
111 | * %$A/g^{m_0} = (g^2)^{m'}$% and we can proceed inductively. The end | |
112 | * result will me %$A/g^m$%. | |
113 | * | |
114 | * Note that this loop keeps track of (what will be) %$r c^m$%, since this | |
115 | * is the result we want, and computes $A/g^m = r^2/a$% on demand. | |
116 | */ | |
9f11b970 | 117 | |
5b295848 MW |
118 | A = mp_sqr(t, r); A = mpmont_reduce(&mm, A, A); |
119 | A = mpmont_mul(&mm, A, A, ainv); /* %$x^t/g^m$% */ | |
9f11b970 | 120 | |
5b295848 MW |
121 | while (s-- > 1) { |
122 | aa = MP_COPY(A); | |
123 | for (i = 1; i < s; i++) | |
124 | { aa = mp_sqr(aa, aa); aa = mpmont_reduce(&mm, aa, aa); } | |
125 | if (!MP_EQ(aa, mm.r)) { | |
9f11b970 | 126 | r = mpmont_mul(&mm, r, r, c); |
5b295848 MW |
127 | A = mp_sqr(A, r); A = mpmont_reduce(&mm, A, A); |
128 | A = mpmont_mul(&mm, A, A, ainv); /* %$x^t/g^m$% */ | |
129 | } | |
130 | c = mp_sqr(c, c); c = mpmont_reduce(&mm, c, c); | |
131 | MP_DROP(aa); | |
9f11b970 | 132 | } |
133 | ||
5b295848 | 134 | /* --- We want the smaller square root --- */ |
9f11b970 | 135 | |
136 | d = mpmont_reduce(&mm, d, r); | |
222c8a43 MW |
137 | r = mp_sub(r, p, d); |
138 | if (MP_CMP(r, <, d)) { mp *tt = r; r = d; d = tt; } | |
5b295848 MW |
139 | |
140 | /* --- Clear away all the temporaries --- */ | |
141 | ||
9f11b970 | 142 | mp_drop(ainv); |
143 | mp_drop(r); mp_drop(c); | |
5b295848 | 144 | mp_drop(A); |
9f11b970 | 145 | mpmont_destroy(&mm); |
146 | ||
5b295848 MW |
147 | /* --- Done --- */ |
148 | ||
9f11b970 | 149 | return (d); |
150 | } | |
151 | ||
152 | /*----- Test rig ----------------------------------------------------------*/ | |
153 | ||
154 | #ifdef TEST_RIG | |
155 | ||
156 | #include <mLib/testrig.h> | |
157 | ||
158 | static int verify(dstr *v) | |
159 | { | |
160 | mp *a = *(mp **)v[0].buf; | |
161 | mp *p = *(mp **)v[1].buf; | |
162 | mp *rr = *(mp **)v[2].buf; | |
163 | mp *r = mp_modsqrt(MP_NEW, a, p); | |
164 | int ok = 0; | |
165 | ||
166 | if (!r) | |
167 | ok = 0; | |
4b536f42 | 168 | else if (MP_EQ(r, rr)) |
9f11b970 | 169 | ok = 1; |
9f11b970 | 170 | |
171 | if (!ok) { | |
172 | fputs("\n*** fail\n", stderr); | |
173 | fputs("a = ", stderr); mp_writefile(a, stderr, 10); fputc('\n', stderr); | |
174 | fputs("p = ", stderr); mp_writefile(p, stderr, 10); fputc('\n', stderr); | |
175 | if (r) { | |
45c0fd36 | 176 | fputs("r = ", stderr); |
9f11b970 | 177 | mp_writefile(r, stderr, 10); |
178 | fputc('\n', stderr); | |
179 | } else | |
45c0fd36 | 180 | fputs("r = <undef>\n", stderr); |
9f11b970 | 181 | fputs("rr = ", stderr); mp_writefile(rr, stderr, 10); fputc('\n', stderr); |
182 | ok = 0; | |
183 | } | |
184 | ||
185 | mp_drop(a); | |
186 | mp_drop(p); | |
f1140c41 | 187 | mp_drop(r); |
9f11b970 | 188 | mp_drop(rr); |
189 | assert(mparena_count(MPARENA_GLOBAL) == 0); | |
190 | return (ok); | |
191 | } | |
192 | ||
193 | static test_chunk tests[] = { | |
194 | { "modsqrt", verify, { &type_mp, &type_mp, &type_mp, 0 } }, | |
195 | { 0, 0, { 0 } } | |
196 | }; | |
197 | ||
198 | int main(int argc, char *argv[]) | |
199 | { | |
200 | sub_init(); | |
0f00dc4c | 201 | test_run(argc, argv, tests, SRCDIR "/t/mp"); |
9f11b970 | 202 | return (0); |
203 | } | |
204 | ||
205 | #endif | |
206 | ||
207 | /*----- That's all, folks -------------------------------------------------*/ |