3 * $Id: mp-modsqrt.c,v 1.4 2001/06/16 12:56:38 mdw Exp $
5 * Compute square roots modulo a prime
7 * (c) 2000 Straylight/Edgeware
10 /*----- Licensing notice --------------------------------------------------*
12 * This file is part of Catacomb.
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
17 * License, or (at your option) any later version.
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.
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,
30 /*----- Revision history --------------------------------------------------*
32 * $Log: mp-modsqrt.c,v $
33 * Revision 1.4 2001/06/16 12:56:38 mdw
34 * Fixes for interface change to @mpmont_expr@ and @mpmont_mexpr@.
36 * Revision 1.3 2001/02/03 12:00:29 mdw
37 * Now @mp_drop@ checks its argument is non-NULL before attempting to free
38 * it. Note that the macro version @MP_DROP@ doesn't do this.
40 * Revision 1.2 2000/10/08 12:02:21 mdw
41 * Use @MP_EQ@ instead of @MP_CMP@.
43 * Revision 1.1 2000/06/22 19:01:31 mdw
44 * Compute square roots in a prime field.
48 /*----- Header files ------------------------------------------------------*/
56 /*----- Main code ---------------------------------------------------------*/
58 /* --- @mp_modsqrt@ --- *
60 * Arguments: @mp *d@ = destination integer
61 * @mp *a@ = source integer
62 * @mp *p@ = modulus (must be prime)
64 * Returns: If %$a$% is a quadratic residue, a square root of %$a$%; else
67 * Use: Returns an integer %$x$% such that %$x^2 \equiv a \pmod{p}$%,
68 * if one exists; else a null pointer. This function will not
69 * work if %$p$% is composite: you must factor the modulus, take
70 * a square root mod each factor, and recombine the results
71 * using the Chinese Remainder Theorem.
74 mp
*mp_modsqrt(mp
*d
, mp
*a
, mp
*p
)
85 /* --- Cope if %$a \not\in Q_p$% --- */
87 if (mp_jacobi(a
, p
) != 1) {
92 /* --- Choose some quadratic non-residue --- */
95 grand
*g
= fibrand_create(0);
99 b
= mprand_range(b
, p
, g
, 0);
100 while (mp_jacobi(b
, p
) != -1);
104 /* --- Find the inverse of %$a$% --- */
107 mp_gcd(0, &ainv
, 0, a
, p
);
109 /* --- Split %$p - 1$% into a power of two and an odd number --- */
111 t
= mp_sub(MP_NEW
, p
, MP_ONE
);
112 t
= mp_odd(t
, t
, &s
);
114 /* --- Now to really get going --- */
116 mpmont_create(&mm
, p
);
117 b
= mpmont_mul(&mm
, b
, b
, mm
.r2
);
118 c
= mpmont_expr(&mm
, b
, b
, t
);
119 t
= mp_add(t
, t
, MP_ONE
);
121 dd
= mpmont_mul(&mm
, MP_NEW
, a
, mm
.r2
);
122 r
= mpmont_expr(&mm
, t
, dd
, t
);
124 ainv
= mpmont_mul(&mm
, ainv
, ainv
, mm
.r2
);
126 mone
= mp_sub(MP_NEW
, p
, mm
.r
);
130 for (i
= 1; i
< s
; i
++) {
132 /* --- Compute %$d_0 = r^2a^{-1}$% --- */
135 dd
= mpmont_reduce(&mm
, dd
, dd
);
136 dd
= mpmont_mul(&mm
, dd
, dd
, ainv
);
138 /* --- Now %$d = d_0^{s - i - 1}$% --- */
140 for (j
= i
; j
< s
- 1; j
++) {
142 dd
= mpmont_reduce(&mm
, dd
, dd
);
145 /* --- Fiddle at the end --- */
148 r
= mpmont_mul(&mm
, r
, r
, c
);
150 c
= mpmont_reduce(&mm
, c
, c
);
153 /* --- Done, so tidy up --- */
155 d
= mpmont_reduce(&mm
, d
, r
);
157 mp_drop(r
); mp_drop(c
);
165 /*----- Test rig ----------------------------------------------------------*/
169 #include <mLib/testrig.h>
171 static int verify(dstr
*v
)
173 mp
*a
= *(mp
**)v
[0].buf
;
174 mp
*p
= *(mp
**)v
[1].buf
;
175 mp
*rr
= *(mp
**)v
[2].buf
;
176 mp
*r
= mp_modsqrt(MP_NEW
, a
, p
);
181 else if (MP_EQ(r
, rr
))
190 fputs("\n*** fail\n", stderr
);
191 fputs("a = ", stderr
); mp_writefile(a
, stderr
, 10); fputc('\n', stderr
);
192 fputs("p = ", stderr
); mp_writefile(p
, stderr
, 10); fputc('\n', stderr
);
194 fputs("r = ", stderr
);
195 mp_writefile(r
, stderr
, 10);
198 fputs("r = <undef>\n", stderr
);
199 fputs("rr = ", stderr
); mp_writefile(rr
, stderr
, 10); fputc('\n', stderr
);
207 assert(mparena_count(MPARENA_GLOBAL
) == 0);
211 static test_chunk tests
[] = {
212 { "modsqrt", verify
, { &type_mp
, &type_mp
, &type_mp
, 0 } },
216 int main(int argc
, char *argv
[])
219 test_run(argc
, argv
, tests
, SRCDIR
"/tests/mp");
225 /*----- That's all, folks -------------------------------------------------*/