- /* --- 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.
+ */