chiark / gitweb /
math/ptab.in: Include the correct Oakley 2048 group!
[catacomb] / math / mp-modsqrt.c
CommitLineData
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
57mp *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
158static 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
193static test_chunk tests[] = {
194 { "modsqrt", verify, { &type_mp, &type_mp, &type_mp, 0 } },
195 { 0, 0, { 0 } }
196};
197
198int 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 -------------------------------------------------*/