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$% --- */
75 if (d
!= a
) mp_drop(d
);
80 /* --- Choose some quadratic non-residue --- */
82 gr
= fibrand_create(0);
84 do b
= mprand_range(b
, p
, gr
, 0); while (mp_jacobi(b
, p
) != -1);
87 /* --- Some initial setup --- */
89 mpmont_create(&mm
, p
);
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$% */
94 b
= mpmont_mul(&mm
, b
, b
, mm
.r2
);
95 c
= mpmont_expr(&mm
, b
, b
, t
); /* %$b^t \bmod p$% */
96 t
= mp_add(t
, t
, MP_ONE
);
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$% */
101 /* --- Now for the main loop --- *
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
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$%.
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.
118 A
= mp_sqr(t
, r
); A
= mpmont_reduce(&mm
, A
, A
);
119 A
= mpmont_mul(&mm
, A
, A
, ainv
); /* %$x^t/g^m$% */
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
)) {
126 r
= mpmont_mul(&mm
, r
, r
, c
);
127 A
= mp_sqr(A
, r
); A
= mpmont_reduce(&mm
, A
, A
);
128 A
= mpmont_mul(&mm
, A
, A
, ainv
); /* %$x^t/g^m$% */
130 c
= mp_sqr(c
, c
); c
= mpmont_reduce(&mm
, c
, c
);
134 /* --- We want the smaller square root --- */
136 d
= mpmont_reduce(&mm
, d
, r
);
138 if (MP_CMP(r
, <, d
)) { mp
*tt
= r
; r
= d
; d
= tt
; }
140 /* --- Clear away all the temporaries --- */
143 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 -------------------------------------------------*/