3 * Compute square roots modulo a prime
5 * (c) 2000 Straylight/Edgeware
8 /*----- Licensing notice --------------------------------------------------*
10 * This file is part of Catacomb.
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.
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.
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,
28 /*----- Header files ------------------------------------------------------*/
36 /*----- Main code ---------------------------------------------------------*/
38 /* --- @mp_modsqrt@ --- *
40 * Arguments: @mp *d@ = destination integer
41 * @mp *a@ = source integer
42 * @mp *p@ = modulus (must be prime)
44 * Returns: If %$a$% is a quadratic residue, a square root of %$a$%; else
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.
53 * We guarantee that the square root returned is the smallest
54 * one (i.e., the `positive' square root).
57 mp
*mp_modsqrt(mp
*d
, mp
*a
, mp
*p
)
68 /* --- Cope if %$a \not\in Q_p$% --- */
70 if (mp_jacobi(a
, p
) != 1) {
75 /* --- Choose some quadratic non-residue --- */
78 grand
*g
= fibrand_create(0);
82 b
= mprand_range(b
, p
, g
, 0);
83 while (mp_jacobi(b
, p
) != -1);
87 /* --- Find the inverse of %$a$% --- */
89 ainv
= mp_modinv(MP_NEW
, a
, p
);
91 /* --- Split %$p - 1$% into a power of two and an odd number --- */
93 t
= mp_sub(MP_NEW
, p
, MP_ONE
);
96 /* --- Now to really get going --- */
98 mpmont_create(&mm
, p
);
99 b
= mpmont_mul(&mm
, b
, b
, mm
.r2
);
100 c
= mpmont_expr(&mm
, b
, b
, t
);
101 t
= mp_add(t
, t
, MP_ONE
);
103 dd
= mpmont_mul(&mm
, MP_NEW
, a
, mm
.r2
);
104 r
= mpmont_expr(&mm
, t
, dd
, t
);
106 ainv
= mpmont_mul(&mm
, ainv
, ainv
, mm
.r2
);
108 mone
= mp_sub(MP_NEW
, p
, mm
.r
);
112 for (i
= 1; i
< s
; i
++) {
114 /* --- Compute %$d_0 = r^2a^{-1}$% --- */
117 dd
= mpmont_reduce(&mm
, dd
, dd
);
118 dd
= mpmont_mul(&mm
, dd
, dd
, ainv
);
120 /* --- Now %$d = d_0^{2^{s - i - 1}}$% --- */
122 for (j
= i
; j
< s
- 1; j
++) {
124 dd
= mpmont_reduce(&mm
, dd
, dd
);
127 /* --- Fiddle at the end --- */
130 r
= mpmont_mul(&mm
, r
, r
, c
);
132 c
= mpmont_reduce(&mm
, c
, c
);
135 /* --- Done, so tidy up --- *
137 * Canonify the answer.
140 d
= mpmont_reduce(&mm
, d
, r
);
142 if (MP_CMP(r
, <, d
)) { mp
*tt
= r
; r
= d
; d
= tt
; }
144 mp_drop(r
); mp_drop(c
);
152 /*----- Test rig ----------------------------------------------------------*/
156 #include <mLib/testrig.h>
158 static int verify(dstr
*v
)
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
);
168 else if (MP_EQ(r
, rr
))
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
);
176 fputs("r = ", stderr
);
177 mp_writefile(r
, stderr
, 10);
180 fputs("r = <undef>\n", stderr
);
181 fputs("rr = ", stderr
); mp_writefile(rr
, stderr
, 10); fputc('\n', stderr
);
189 assert(mparena_count(MPARENA_GLOBAL
) == 0);
193 static test_chunk tests
[] = {
194 { "modsqrt", verify
, { &type_mp
, &type_mp
, &type_mp
, 0 } },
198 int main(int argc
, char *argv
[])
201 test_run(argc
, argv
, tests
, SRCDIR
"/t/mp");
207 /*----- That's all, folks -------------------------------------------------*/