1 /* -*-c-*-
2 *
3 * $Id: mp-modsqrt.c,v 1.5 2004/04/08 01:36:15 mdw Exp$
4 *
5 * Compute square roots modulo a prime
6 *
7 * (c) 2000 Straylight/Edgeware
8 */
10 /*----- Licensing notice --------------------------------------------------*
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
18 *
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.
23 *
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 */
32 #include "fibrand.h"
33 #include "grand.h"
34 #include "mp.h"
35 #include "mpmont.h"
36 #include "mprand.h"
38 /*----- Main code ---------------------------------------------------------*/
40 /* --- @mp_modsqrt@ --- *
41 *
42 * Arguments: @mp *d@ = destination integer
43 * @mp *a@ = source integer
44 * @mp *p@ = modulus (must be prime)
45 *
46 * Returns: If %$a$% is a quadratic residue, a square root of %$a$%; else
47 * a null pointer.
48 *
49 * Use: Returns an integer %$x$% such that %$x^2 \equiv a \pmod{p}$%,
50 * if one exists; else a null pointer. This function will not
51 * work if %$p$% is composite: you must factor the modulus, take
52 * a square root mod each factor, and recombine the results
53 * using the Chinese Remainder Theorem.
54 */
56 mp *mp_modsqrt(mp *d, mp *a, mp *p)
57 {
58 mpmont mm;
59 mp *t;
60 size_t s;
61 mp *b;
62 mp *ainv;
63 mp *c, *r;
64 size_t i, j;
65 mp *dd, *mone;
67 /* --- Cope if %$a \not\in Q_p$% --- */
69 if (mp_jacobi(a, p) != 1) {
70 mp_drop(d);
71 return (0);
72 }
74 /* --- Choose some quadratic non-residue --- */
76 {
77 grand *g = fibrand_create(0);
79 b = MP_NEW;
80 do
81 b = mprand_range(b, p, g, 0);
82 while (mp_jacobi(b, p) != -1);
83 g->ops->destroy(g);
84 }
86 /* --- Find the inverse of %$a$% --- */
88 ainv = mp_modinv(MP_NEW, a, p);
90 /* --- Split %$p - 1$% into a power of two and an odd number --- */
92 t = mp_sub(MP_NEW, p, MP_ONE);
93 t = mp_odd(t, t, &s);
95 /* --- Now to really get going --- */
97 mpmont_create(&mm, p);
98 b = mpmont_mul(&mm, b, b, mm.r2);
99 c = mpmont_expr(&mm, b, b, t);
100 t = mp_add(t, t, MP_ONE);
101 t = mp_lsr(t, t, 1);
102 dd = mpmont_mul(&mm, MP_NEW, a, mm.r2);
103 r = mpmont_expr(&mm, t, dd, t);
104 mp_drop(dd);
105 ainv = mpmont_mul(&mm, ainv, ainv, mm.r2);
107 mone = mp_sub(MP_NEW, p, mm.r);
109 dd = MP_NEW;
111 for (i = 1; i < s; i++) {
113 /* --- Compute %$d_0 = r^2a^{-1}$% --- */
115 dd = mp_sqr(dd, r);
116 dd = mpmont_reduce(&mm, dd, dd);
117 dd = mpmont_mul(&mm, dd, dd, ainv);
119 /* --- Now %$d = d_0^{s - i - 1}$% --- */
121 for (j = i; j < s - 1; j++) {
122 dd = mp_sqr(dd, dd);
123 dd = mpmont_reduce(&mm, dd, dd);
124 }
126 /* --- Fiddle at the end --- */
128 if (MP_EQ(dd, mone))
129 r = mpmont_mul(&mm, r, r, c);
130 c = mp_sqr(c, c);
131 c = mpmont_reduce(&mm, c, c);
132 }
134 /* --- Done, so tidy up --- */
136 d = mpmont_reduce(&mm, d, r);
137 mp_drop(ainv);
138 mp_drop(r); mp_drop(c);
139 mp_drop(dd);
140 mp_drop(mone);
141 mpmont_destroy(&mm);
143 return (d);
144 }
146 /*----- Test rig ----------------------------------------------------------*/
148 #ifdef TEST_RIG
150 #include <mLib/testrig.h>
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;
160 if (!r)
161 ok = 0;
162 else if (MP_EQ(r, rr))
163 ok = 1;
164 else {
165 r = mp_sub(r, p, r);
166 if (MP_EQ(r, rr))
167 ok = 1;
168 }
170 if (!ok) {
171 fputs("\n*** fail\n", stderr);
172 fputs("a = ", stderr); mp_writefile(a, stderr, 10); fputc('\n', stderr);
173 fputs("p = ", stderr); mp_writefile(p, stderr, 10); fputc('\n', stderr);
174 if (r) {
175 fputs("r = ", stderr);
176 mp_writefile(r, stderr, 10);
177 fputc('\n', stderr);
178 } else
179 fputs("r = <undef>\n", stderr);
180 fputs("rr = ", stderr); mp_writefile(rr, stderr, 10); fputc('\n', stderr);
181 ok = 0;
182 }
184 mp_drop(a);
185 mp_drop(p);
186 mp_drop(r);
187 mp_drop(rr);
188 assert(mparena_count(MPARENA_GLOBAL) == 0);
189 return (ok);
190 }
192 static test_chunk tests[] = {
193 { "modsqrt", verify, { &type_mp, &type_mp, &type_mp, 0 } },
194 { 0, 0, { 0 } }
195 };
197 int main(int argc, char *argv[])
198 {
199 sub_init();
200 test_run(argc, argv, tests, SRCDIR "/tests/mp");
201 return (0);
202 }
204 #endif
206 /*----- That's all, folks -------------------------------------------------*/