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
)
67 /* --- Cope if %$a \not\in Q_p$% --- */
69 if (mp_jacobi(a
, p
) != 1) {
74 /* --- Choose some quadratic non-residue --- */
76 gr
= fibrand_create(0);
78 do b
= mprand_range(b
, p
, gr
, 0); while (mp_jacobi(b
, p
) != -1);
81 /* --- Some initial setup --- */
83 mpmont_create(&mm
, p
);
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$% */
88 b
= mpmont_mul(&mm
, b
, b
, mm
.r2
);
89 c
= mpmont_expr(&mm
, b
, b
, t
); /* %$b^t \bmod p$% */
90 t
= mp_add(t
, t
, MP_ONE
);
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$% */
95 /* --- Now for the main loop --- *
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
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$%.
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.
112 A
= mp_sqr(t
, r
); A
= mpmont_reduce(&mm
, A
, A
);
113 A
= mpmont_mul(&mm
, A
, A
, ainv
); /* %$x^t/g^m$% */
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
)) {
120 r
= mpmont_mul(&mm
, r
, r
, c
);
121 A
= mp_sqr(A
, r
); A
= mpmont_reduce(&mm
, A
, A
);
122 A
= mpmont_mul(&mm
, A
, A
, ainv
); /* %$x^t/g^m$% */
124 c
= mp_sqr(c
, c
); c
= mpmont_reduce(&mm
, c
, c
);
128 /* --- We want the smaller square root --- */
130 d
= mpmont_reduce(&mm
, d
, r
);
132 if (MP_CMP(r
, <, d
)) { mp
*tt
= r
; r
= d
; d
= tt
; }
134 /* --- Clear away all the temporaries --- */
137 mp_drop(r
); mp_drop(c
);
146 /*----- Test rig ----------------------------------------------------------*/
150 #include <mLib/testrig.h>
152 static int verify(dstr
*v
)
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
);
162 else if (MP_EQ(r
, rr
))
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
);
170 fputs("r = ", stderr
);
171 mp_writefile(r
, stderr
, 10);
174 fputs("r = <undef>\n", stderr
);
175 fputs("rr = ", stderr
); mp_writefile(rr
, stderr
, 10); fputc('\n', stderr
);
183 assert(mparena_count(MPARENA_GLOBAL
) == 0);
187 static test_chunk tests
[] = {
188 { "modsqrt", verify
, { &type_mp
, &type_mp
, &type_mp
, 0 } },
192 int main(int argc
, char *argv
[])
195 test_run(argc
, argv
, tests
, SRCDIR
"/t/mp");
201 /*----- That's all, folks -------------------------------------------------*/