mp *mp_modsqrt(mp *d, mp *a, mp *p)
{
mpmont mm;
- mp *t;
- size_t s;
- mp *b;
+ size_t i, s;
+ mp *b, *c;
mp *ainv;
- mp *c, *r;
- size_t i, j;
- mp *dd, *mone;
+ mp *r, *A, *aa;
+ mp *t;
+ grand *gr;
+ int j;
/* --- Cope if %$a \not\in Q_p$% --- */
- if (mp_jacobi(a, p) != 1) {
+ j = mp_jacobi(a, p);
+ if (j == -1) {
mp_drop(d);
return (0);
+ } else if (j == 0) {
+ if (d != a) mp_drop(d);
+ d = MP_COPY(a);
+ return (d);
}
/* --- Choose some quadratic non-residue --- */
- {
- grand *g = fibrand_create(0);
-
- b = MP_NEW;
- do
- b = mprand_range(b, p, g, 0);
- while (mp_jacobi(b, p) != -1);
- g->ops->destroy(g);
- }
-
- /* --- Find the inverse of %$a$% --- */
-
- ainv = mp_modinv(MP_NEW, a, p);
+ gr = fibrand_create(0);
+ b = MP_NEW;
+ do b = mprand_range(b, p, gr, 0); while (mp_jacobi(b, p) != -1);
+ gr->ops->destroy(gr);
- /* --- Split %$p - 1$% into a power of two and an odd number --- */
-
- t = mp_sub(MP_NEW, p, MP_ONE);
- t = mp_odd(t, t, &s);
-
- /* --- Now to really get going --- */
+ /* --- Some initial setup --- */
mpmont_create(&mm, p);
+ ainv = mp_modinv(MP_NEW, a, p); /* %$a^{-1} \bmod p$% */
+ ainv = mpmont_mul(&mm, ainv, ainv, mm.r2);
+ t = mp_sub(MP_NEW, p, MP_ONE);
+ t = mp_odd(t, t, &s); /* %$2^s t = p - 1$% */
b = mpmont_mul(&mm, b, b, mm.r2);
- c = mpmont_expr(&mm, b, b, t);
+ c = mpmont_expr(&mm, b, b, t); /* %$b^t \bmod p$% */
t = mp_add(t, t, MP_ONE);
- t = mp_lsr(t, t, 1);
- dd = mpmont_mul(&mm, MP_NEW, a, mm.r2);
- r = mpmont_expr(&mm, t, dd, t);
- mp_drop(dd);
- ainv = mpmont_mul(&mm, ainv, ainv, mm.r2);
-
- mone = mp_sub(MP_NEW, p, mm.r);
-
- dd = MP_NEW;
-
- for (i = 1; i < s; i++) {
+ t = mp_lsr(t, t, 1); /* %$(t + 1)/2$% */
+ a = mpmont_mul(&mm, MP_NEW, a, mm.r2);
+ r = mpmont_expr(&mm, a, a, t); /* %$a^{(t+1)/2} \bmod p$% */
- /* --- Compute %$d_0 = r^2a^{-1}$% --- */
-
- dd = mp_sqr(dd, r);
- dd = mpmont_reduce(&mm, dd, dd);
- dd = mpmont_mul(&mm, dd, dd, ainv);
-
- /* --- Now %$d = d_0^{2^{s - i - 1}}$% --- */
-
- for (j = i; j < s - 1; j++) {
- dd = mp_sqr(dd, dd);
- dd = mpmont_reduce(&mm, dd, dd);
- }
+ /* --- Now for the main loop --- *
+ *
+ * Let %$g = c^{-2}$%; we know that %$g$% is a generator of the order-
+ * %$2^{s-1}$% subgroup mod %$p$%. We also know that %$A = a^t = r^2/a$%
+ * is an element of this group. If we can determine %$m$% such that
+ * %$g^m = A$% then %$a^{(t+1)/2}/g^{m/2} = r c^m$% is the square root we
+ * seek.
+ *
+ * Write %$m = m_0 + 2 m'$%. Then %$A^{2^{s-1}} = g^{m_0 2^{s-1}}$%, which
+ * is %$1$% if %$m_0 = 0$% or %$-1$% if %$m_0 = 1$% (modulo %$p$%). Then
+ * %$A/g^{m_0} = (g^2)^{m'}$% and we can proceed inductively. The end
+ * result will me %$A/g^m$%.
+ *
+ * Note that this loop keeps track of (what will be) %$r c^m$%, since this
+ * is the result we want, and computes $A/g^m = r^2/a$% on demand.
+ */
- /* --- Fiddle at the end --- */
+ A = mp_sqr(t, r); A = mpmont_reduce(&mm, A, A);
+ A = mpmont_mul(&mm, A, A, ainv); /* %$x^t/g^m$% */
- if (MP_EQ(dd, mone))
+ while (s-- > 1) {
+ aa = MP_COPY(A);
+ for (i = 1; i < s; i++)
+ { aa = mp_sqr(aa, aa); aa = mpmont_reduce(&mm, aa, aa); }
+ if (!MP_EQ(aa, mm.r)) {
r = mpmont_mul(&mm, r, r, c);
- c = mp_sqr(c, c);
- c = mpmont_reduce(&mm, c, c);
+ A = mp_sqr(A, r); A = mpmont_reduce(&mm, A, A);
+ A = mpmont_mul(&mm, A, A, ainv); /* %$x^t/g^m$% */
+ }
+ c = mp_sqr(c, c); c = mpmont_reduce(&mm, c, c);
+ MP_DROP(aa);
}
- /* --- Done, so tidy up --- *
- *
- * Canonify the answer.
- */
+ /* --- We want the smaller square root --- */
d = mpmont_reduce(&mm, d, r);
r = mp_sub(r, p, d);
if (MP_CMP(r, <, d)) { mp *tt = r; r = d; d = tt; }
+
+ /* --- Clear away all the temporaries --- */
+
mp_drop(ainv);
mp_drop(r); mp_drop(c);
- mp_drop(dd);
- mp_drop(mone);
+ mp_drop(A);
mpmont_destroy(&mm);
+ /* --- Done --- */
+
return (d);
}