Commit | Line | Data |
---|---|---|
9f11b970 | 1 | /* -*-c-*- |
2 | * | |
9f11b970 | 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; | |
9f11b970 | 66 | |
67 | /* --- Cope if %$a \not\in Q_p$% --- */ | |
68 | ||
69 | if (mp_jacobi(a, p) != 1) { | |
f1140c41 | 70 | mp_drop(d); |
9f11b970 | 71 | return (0); |
72 | } | |
73 | ||
74 | /* --- Choose some quadratic non-residue --- */ | |
75 | ||
5b295848 MW |
76 | gr = fibrand_create(0); |
77 | b = MP_NEW; | |
78 | do b = mprand_range(b, p, gr, 0); while (mp_jacobi(b, p) != -1); | |
79 | gr->ops->destroy(gr); | |
45c0fd36 | 80 | |
5b295848 | 81 | /* --- Some initial setup --- */ |
9f11b970 | 82 | |
83 | mpmont_create(&mm, p); | |
5b295848 MW |
84 | ainv = mp_modinv(MP_NEW, a, p); /* %$a^{-1} \bmod p$% */ |
85 | ainv = mpmont_mul(&mm, ainv, ainv, mm.r2); | |
86 | t = mp_sub(MP_NEW, p, MP_ONE); | |
87 | t = mp_odd(t, t, &s); /* %$2^s t = p - 1$% */ | |
b0b682aa | 88 | b = mpmont_mul(&mm, b, b, mm.r2); |
5b295848 | 89 | c = mpmont_expr(&mm, b, b, t); /* %$b^t \bmod p$% */ |
9f11b970 | 90 | t = mp_add(t, t, MP_ONE); |
5b295848 MW |
91 | t = mp_lsr(t, t, 1); /* %$(t + 1)/2$% */ |
92 | a = mpmont_mul(&mm, MP_NEW, a, mm.r2); | |
93 | r = mpmont_expr(&mm, a, a, t); /* %$a^{(t+1)/2} \bmod p$% */ | |
9f11b970 | 94 | |
5b295848 MW |
95 | /* --- Now for the main loop --- * |
96 | * | |
97 | * Let %$g = c^{-2}$%; we know that %$g$% is a generator of the order- | |
98 | * %$2^{s-1}$% subgroup mod %$p$%. We also know that %$A = a^t = r^2/a$% | |
99 | * is an element of this group. If we can determine %$m$% such that | |
100 | * %$g^m = A$% then %$a^{(t+1)/2}/g^{m/2} = r c^m$% is the square root we | |
101 | * seek. | |
102 | * | |
103 | * Write %$m = m_0 + 2 m'$%. Then %$A^{2^{s-1}} = g^{m_0 2^{s-1}}$%, which | |
104 | * is %$1$% if %$m_0 = 0$% or %$-1$% if %$m_0 = 1$% (modulo %$p$%). Then | |
105 | * %$A/g^{m_0} = (g^2)^{m'}$% and we can proceed inductively. The end | |
106 | * result will me %$A/g^m$%. | |
107 | * | |
108 | * Note that this loop keeps track of (what will be) %$r c^m$%, since this | |
109 | * is the result we want, and computes $A/g^m = r^2/a$% on demand. | |
110 | */ | |
9f11b970 | 111 | |
5b295848 MW |
112 | A = mp_sqr(t, r); A = mpmont_reduce(&mm, A, A); |
113 | A = mpmont_mul(&mm, A, A, ainv); /* %$x^t/g^m$% */ | |
9f11b970 | 114 | |
5b295848 MW |
115 | while (s-- > 1) { |
116 | aa = MP_COPY(A); | |
117 | for (i = 1; i < s; i++) | |
118 | { aa = mp_sqr(aa, aa); aa = mpmont_reduce(&mm, aa, aa); } | |
119 | if (!MP_EQ(aa, mm.r)) { | |
9f11b970 | 120 | r = mpmont_mul(&mm, r, r, c); |
5b295848 MW |
121 | A = mp_sqr(A, r); A = mpmont_reduce(&mm, A, A); |
122 | A = mpmont_mul(&mm, A, A, ainv); /* %$x^t/g^m$% */ | |
123 | } | |
124 | c = mp_sqr(c, c); c = mpmont_reduce(&mm, c, c); | |
125 | MP_DROP(aa); | |
9f11b970 | 126 | } |
127 | ||
5b295848 | 128 | /* --- We want the smaller square root --- */ |
9f11b970 | 129 | |
130 | d = mpmont_reduce(&mm, d, r); | |
222c8a43 MW |
131 | r = mp_sub(r, p, d); |
132 | if (MP_CMP(r, <, d)) { mp *tt = r; r = d; d = tt; } | |
5b295848 MW |
133 | |
134 | /* --- Clear away all the temporaries --- */ | |
135 | ||
9f11b970 | 136 | mp_drop(ainv); |
137 | mp_drop(r); mp_drop(c); | |
5b295848 | 138 | mp_drop(A); |
9f11b970 | 139 | mpmont_destroy(&mm); |
140 | ||
5b295848 MW |
141 | /* --- Done --- */ |
142 | ||
9f11b970 | 143 | return (d); |
144 | } | |
145 | ||
146 | /*----- Test rig ----------------------------------------------------------*/ | |
147 | ||
148 | #ifdef TEST_RIG | |
149 | ||
150 | #include <mLib/testrig.h> | |
151 | ||
152 | static int verify(dstr *v) | |
153 | { | |
154 | mp *a = *(mp **)v[0].buf; | |
155 | mp *p = *(mp **)v[1].buf; | |
156 | mp *rr = *(mp **)v[2].buf; | |
157 | mp *r = mp_modsqrt(MP_NEW, a, p); | |
158 | int ok = 0; | |
159 | ||
160 | if (!r) | |
161 | ok = 0; | |
4b536f42 | 162 | else if (MP_EQ(r, rr)) |
9f11b970 | 163 | ok = 1; |
9f11b970 | 164 | |
165 | if (!ok) { | |
166 | fputs("\n*** fail\n", stderr); | |
167 | fputs("a = ", stderr); mp_writefile(a, stderr, 10); fputc('\n', stderr); | |
168 | fputs("p = ", stderr); mp_writefile(p, stderr, 10); fputc('\n', stderr); | |
169 | if (r) { | |
45c0fd36 | 170 | fputs("r = ", stderr); |
9f11b970 | 171 | mp_writefile(r, stderr, 10); |
172 | fputc('\n', stderr); | |
173 | } else | |
45c0fd36 | 174 | fputs("r = <undef>\n", stderr); |
9f11b970 | 175 | fputs("rr = ", stderr); mp_writefile(rr, stderr, 10); fputc('\n', stderr); |
176 | ok = 0; | |
177 | } | |
178 | ||
179 | mp_drop(a); | |
180 | mp_drop(p); | |
f1140c41 | 181 | mp_drop(r); |
9f11b970 | 182 | mp_drop(rr); |
183 | assert(mparena_count(MPARENA_GLOBAL) == 0); | |
184 | return (ok); | |
185 | } | |
186 | ||
187 | static test_chunk tests[] = { | |
188 | { "modsqrt", verify, { &type_mp, &type_mp, &type_mp, 0 } }, | |
189 | { 0, 0, { 0 } } | |
190 | }; | |
191 | ||
192 | int main(int argc, char *argv[]) | |
193 | { | |
194 | sub_init(); | |
0f00dc4c | 195 | test_run(argc, argv, tests, SRCDIR "/t/mp"); |
9f11b970 | 196 | return (0); |
197 | } | |
198 | ||
199 | #endif | |
200 | ||
201 | /*----- That's all, folks -------------------------------------------------*/ |