From 5b295848664fa50612fb2df8aba58f29f9f74d63 Mon Sep 17 00:00:00 2001
From: Mark Wooding <mdw@distorted.org.uk>
Date: Sat, 3 Aug 2013 12:28:45 +0100
Subject: [PATCH] math/mp-modsqrt.c: Reformat and add commentary.

---
 math/mp-modsqrt.c | 112 ++++++++++++++++++++++++++----------------------------
 1 file changed, 53 insertions(+), 59 deletions(-)

diff --git a/math/mp-modsqrt.c b/math/mp-modsqrt.c
index b4a11fa..c5aed1f 100644
--- a/math/mp-modsqrt.c
+++ b/math/mp-modsqrt.c
@@ -57,13 +57,12 @@
 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;
 
   /* --- Cope if %$a \not\in Q_p$% --- */
 
@@ -74,78 +73,73 @@ mp *mp_modsqrt(mp *d, mp *a, mp *p)
 
   /* --- 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);
 }
 
-- 
2.11.0