From fa17e5dc48bf2961e5a99c83ed3d8018b329d037 Mon Sep 17 00:00:00 2001
From: Mark Wooding <mdw@distorted.org.uk>
Date: Fri, 9 Sep 2016 11:06:41 +0100
Subject: [PATCH] math/{mpbarrett,mpmont}.h: Provide correctness proofs for
 these methods.

Add commentary explaining how these reduction algorithms actually work,
with proofs.
---
 math/mpbarrett.h | 22 ++++++++++++++++++++++
 math/mpmont.h    | 34 ++++++++++++++++++++++++----------
 2 files changed, 46 insertions(+), 10 deletions(-)

diff --git a/math/mpbarrett.h b/math/mpbarrett.h
index 40e0fe42..b3ab798e 100644
--- a/math/mpbarrett.h
+++ b/math/mpbarrett.h
@@ -40,6 +40,28 @@
  * modexp routine provided which uses Barrett reduction rather than
  * Montgomery reduction.  This is handy when you're working on indices in an
  * even-order cyclic group or something.
+ *
+ * In more detail: suppose that %$b^{k-1} \le m < b^k$%.  Let %$\mu = {}$%
+ * %$\lfloor b^{2k}/m \rfloor$%; %$\mu$% is a scaled approximation to the
+ * reciprocal %$1/m$%.  Now, suppose we're given some %$a$% with
+ * %$0 \le a < b^{2k}$%.  The first step is to calculate an approximation
+ * %$q = \lfloor \mu \lfloor a/b^{k-1} \rfloor/b^{k+1} \rfloor$% to the
+ * quotient %$a/m$%.  Then we have:
+ *
+ *	%$\lfloor a/m - a/b^{2k} - b^{k-1}/m + 1/b^{k+1} \rfloor \le {}$%
+ *		%$q \le \lfloor a/m \rfloor
+ *
+ * But by assumption %$a < b^{2k}$% and %$2^{k-1} \le m$% so
+ *
+ *	%$\lfloor a/m \rfloor - 2 \le q \le \lfloor a/m \rfloor$%
+ *
+ * Now we approximate the remainder by calculating %$r = a - q m$%.
+ * Certainly %$r \equiv a \pmod{m}$%; and
+ *
+ *	%$0 \le r \le (a - m \lfloor a/m \rfloor) + 2 m < 3 m$%.
+ *
+ * So the remainder can be fixed up with at most two conditional
+ * subtractions.
  */
 
 #ifndef CATACOMB_MPBARRETT_H
diff --git a/math/mpmont.h b/math/mpmont.h
index 9973ec3f..b0c119af 100644
--- a/math/mpmont.h
+++ b/math/mpmont.h
@@ -54,22 +54,36 @@
  *   * %$b$%, the radix of the number system you're in (here, it's
  *     @MPW_MAX + 1@).
  *
- *   * %$-m^{-1} \bmod b$%, a useful number for the reduction step.  (This
- *     means that the modulus mustn't be even.  This shouldn't be a problem.)
+ *   * %$m' = -m^{-1} \bmod b$%, a useful number for the reduction step.
+ *     (This means that the modulus mustn't be even.  This shouldn't be a
+ *     problem.)
  *
  *   * %$R = b^n > m > b^{n - 1}$%, or at least %$\log_2 R$%.
  *
  *   * %$R \bmod m$% and %$R^2 \bmod m$%, which are useful when doing
  *     calculations such as exponentiation.
  *
- * The result of a Montgomery reduction of %$x$% is %$x R^{-1} \bmod m$%,
- * which doesn't look ever-so useful.  The trick is to initially apply a
- * factor of %$R$% to all of your numbers so that when you multiply and
- * perform a Montgomery reduction you get %$(x R \cdot y R) R^{-1} \bmod m$%,
- * which is just %$x y R \bmod m$%.  Thanks to distributivity, even additions
- * and subtractions can be performed on numbers in this form -- the extra
- * factor of %$R$% just runs through all the calculations until it's finally
- * stripped out by a final reduction operation.
+ * Suppose that %$0 \le a_i \le (b^n + b^i - 1) m$% with %$a_i \equiv {}$%
+ * %$0 \pmod{b^i}$%.  Let %$w_i = m' a_i/b^i \bmod b$%, and set %$a_{i+1} =
+ * a_i + b^i w_i m$%.  Then obviously %$a_{i+1} \equiv {} $% %$a_i
+ * \pmod{m}$%, and less obviously %$a_{i+1}/b^i \equiv a_i/b^i + {}$% %$m m'
+ * a_i/b^i \equiv 0 \pmod{b}$% so %$a_{i+1} \equiv 0 \pmod{b^{i+1}}$%.
+ * Finally, we can bound %$a_{i+1} \le {}$% %$a_i + b^i (b - 1) m = {}$%
+ * %$a_i + (b^{i+1} - b^i) m \le (b^n + b^{i+1} - 1) m$%.  As a result, if
+ * we're given some %a_0%, we can calculate %$a_n \equiv 0 \pmod{R}$%, with
+ * $%a_n \equiv a_0 \pmod{n}$%, i.e., %$a_n/R \equiv a_0 R^{-1} \pmod{m}$%;
+ * furthermore, if %$0 \le a_0 < m + b^n%$ then %$0 \le a_n/R < 2 m$%, so a
+ * fully reduced result can be obtained with a single conditional
+ * subtraction.
+ *
+ * The result of reduing %$a$% is then %$a R^{-1}$% \bmod m$%.  This is
+ * actually rather useful for reducing products, if we run an extra factor of
+ * %$R$% through the calculation: the result of reducing the product of
+ * %$(x R)(y R) = x y R^2$% is then %$x y R \bmod m$%, preserving the running
+ * factor.  Thanks to distributivity, additions and subtractions can be
+ * performed on numbers in this form -- the extra factor of %$R$% just runs
+ * through all the calculations until it's finally stripped out by a final
+ * reduction operation.
  */
 
 /*----- Data structures ---------------------------------------------------*/
-- 
2.11.0