math/: Improve some commentary in the binary-field arithmetic.

  * Explain why `gfreduce_trace' can safely return its answer as an int.

  * Explain how `gfreduce_quadsolve' actually works.  Also explicitly
    guarantee that its result is deterministic.

  * Explain how the `find' method works in `ec-bin.c'.

There's a little fiddling with braces to fit the new commentary in, but
no significant code change.

diff --git a/math/ec-bin.c b/math/ec-bin.c
index d91b034..c7fe96d 100644 (file)
--- a/math/ec-bin.c
+++ b/math/ec-bin.c
@@ -73,6 +73,7 @@ static ec *ecfind(ec_curve *c, ec *d, mp *x)
       v = F_MUL(f, v, u, y);       /* %$B = A x^{-2} = x + a + b x^{-2}$% */
       y = F_QUADSOLVE(f, y, v);                /* %$z^2 + z = B$% */
       if (y) y = F_MUL(f, y, y, x);    /* %$y = z x$% */
+      /* Hence %$y^2 + x y = (z^2 + z) x^2 = A$% */
     }
     MP_DROP(u);
     MP_DROP(v);

diff --git a/math/gfreduce.c b/math/gfreduce.c
index 582b723..83ffeb1 100644 (file)
--- a/math/gfreduce.c
+++ b/math/gfreduce.c
@@ -406,6 +406,12 @@ mp *gfreduce_sqrt(gfreduce *r, mp *d, mp *x)
   unsigned long m = mp_bits(r->p) - 1;
   unsigned long i;
 
+  /* --- This is pretty easy --- *
+   *
+   * Note that %$x = x^{2^m}$%; therefore %$(x^{2^{m-1}})^2 = x^{2^m} = x$%,
+   * so %$x^{2^{m-1}}$% is the square root we seek.
+   */
+
   for (i = 0; i < m - 1; i++) {
     mp *t = gf_sqr(spare, y);
     spare = y;
@@ -428,7 +434,10 @@ mp *gfreduce_sqrt(gfreduce *r, mp *d, mp *x)
  *             @mp *x@ = some polynomial
  *
  * Returns:    The trace of @x@. (%$\Tr(x)=x + x^2 + \cdots + x^{2^{m-1}}$%
- *             if %$x \in \gf{2^m}$%).
+ *             if %$x \in \gf{2^m}$%).  Since the trace is invariant under
+ *             the Frobenius automorphism (i.e., %$\Tr(x)^2 = \Tr(x)$%), it
+ *             must be an element of the base field, i.e., %$\gf{2}$%, and
+ *             we only need a single bit to represent it.
  */
 
 int gfreduce_trace(gfreduce *r, mp *x)
@@ -489,7 +498,18 @@ mp *gfreduce_halftrace(gfreduce *r, mp *d, mp *x)
  *             @mp *d@ = destination
  *             @mp *x@ = some polynomial
  *
- * Returns:    A polynomial @y@ such that %$y^2 + y = x$%, or null.
+ * Returns:    A polynomial @z@ such that %$z^2 + z = x$%, or null.
+ *
+ * Use:                Solves quadratic equations in a field with characteristic 2.
+ *             Suppose we have an equation %$y^2 + A y + B = 0$% where
+ *             %$A \ne 0$%.  (If %$A = 0$% then %$y = \sqrt{B}$% and you
+ *             want @gfreduce_sqrt@ instead.)  Use this function to solve
+ *             %$z^2 + z = B/A^2$%; then set %$y = A z$%, since
+ *             %$y^2 + y = A^2 z^2 + A^2 z = A^2 (z^2 + z) = B$% as
+ *             required.
+ *
+ *             The two roots are %$z$% and %$z + 1$%; this function always
+ *             returns the one with zero scalar coefficient.
  */
 
 mp *gfreduce_quadsolve(gfreduce *r, mp *d, mp *x)
@@ -497,15 +517,55 @@ mp *gfreduce_quadsolve(gfreduce *r, mp *d, mp *x)
   unsigned long m = mp_bits(r->p) - 1;
   mp *t;
 
+  /* --- About the solutions --- *
+   *
+   * Factor %$z^2 + z = z (z + 1)$%.  Therefore, if %$z^2 + z = x$% and
+   * %$z' = z + 1$% then %$z'^2 + z' = z^2 + 1 + z + 1 = z^2 + z$%, so
+   * %$z + 1$% is the other solution.
+   *
+   * A solution exists if and only if %$\Tr(x) = 0$%.  To see the `only if'
+   * implication, recall that the trace function is linear, and hence
+   * $%\Tr(z^2 + z) = \Tr(z)^2 + \Tr(z) = \Tr(z) + \Tr(z) = 0$%.  The `if'
+   * direction will be proven using explicit constructions captured in the
+   * code below.
+   */
+
   MP_COPY(x);
-  if (m & 1)
+  if (m & 1) {
+
+    /* --- A short-cut for fields with odd degree ---
+     *
+     * The method below works in all binary fields, but there's a quicker way
+     * which works whenever the degree is odd.  The half-trace is
+     * %$z = \sum_{0\le i\le (m-1)/2} x^{2^{2i}}$%.  Then %$z^2 + z = {}$%
+     * %$\sum_{0\le i\le (m-1)/2} (x^{2^{2i}} + x^{2^{2i+1}}) = {}$%
+     * %$\Tr(x) + x^{2^m} = \Tr(x) + x$%.  This therefore gives us the
+     * solution we want whenever %$\Tr(x) = 0$%.
+     */
+
     d = gfreduce_halftrace(r, d, x);
-  else {
+  } else {
     mp *z, *w, *rho = MP_NEW;
     mp *spare = MP_NEW;
     grand *fr = fibrand_create(0);
     unsigned long i;
 
+    /* --- Unpicking the magic --- *
+     *
+     * Choose %$\rho \inr \gf{2^m}$% with %$\Tr(\rho) = 1$%.  Let
+     * %$z = \sum_{0\le i<m} \rho^{2^i} \sum_{0\le j<i} x^{2^j} = {}$%
+     * %$\rho^2 x + \rho^4 (x + x^2) + \rho^8 (x + x^2 + x^4) + \cdots + {}$%
+     * %$\rho^{2^{m-1}} (x + x^2 + x^{2^{m-2}})$%.  Then %$z^2 = {}$%
+     * %$\sum_{0\le i<m} \rho^{2^{i+1}} \sum_{0\le j<i} x^{2^{j+1}} = {}$%
+     * %$\sum_{1\le i\le m} \rho^{2^i} \sum_{1\le j\le i} x^{2^j}$% and,
+     * somewhat miraculously, %$z^2 + z = \sum_{0\le i<m} \rho^{2^i} x + {}$%
+     * %$\rho \sum_{1\le i<m} x^{2^i} = x \Tr(\rho) + \rho \Tr(x)$%.  Again,
+     * this gives us the root we want whenever %$\Tr(x) = 0$%.
+     *
+     * The loop below calculates %$w = \Tr(\rho)$% and %$z$% simultaneously,
+     * since the same powers of %$\rho$% are wanted in both calculations.
+     */
+
     for (;;) {
       rho = mprand(rho, m, fr, 0);
       z = MP_ZERO;
@@ -530,6 +590,12 @@ mp *gfreduce_quadsolve(gfreduce *r, mp *d, mp *x)
     d = z;
   }
 
+  /* --- Check that we calculated the right answer --- *
+   *
+   * It should be correct; if it's not then maybe the ring we're working in
+   * isn't really a field.
+   */
+
   t = gf_sqr(MP_NEW, d); t = gfreduce_do(r, t, t); t = gf_add(t, t, d);
   if (!MP_EQ(t, x)) {
     MP_DROP(d);
@@ -537,6 +603,13 @@ mp *gfreduce_quadsolve(gfreduce *r, mp *d, mp *x)
   }
   MP_DROP(t);
   MP_DROP(x);
+
+  /* --- Pick a canonical root --- *
+   *
+   * The two roots are %$z$% and %$z + 1$%; pick the one with a zero
+   * scalar coefficient just for consistency's sake.
+   */
+
   if (d) d->v[0] &= ~(mpw)1;
   return (d);
 }

diff --git a/math/gfreduce.h b/math/gfreduce.h
index 9f69585..ef075c4 100644 (file)
--- a/math/gfreduce.h
+++ b/math/gfreduce.h
@@ -129,7 +129,10 @@ extern mp *gfreduce_sqrt(gfreduce */*r*/, mp */*d*/, mp */*x*/);
  *             @mp *x@ = some polynomial
  *
  * Returns:    The trace of @x@. (%$\Tr(x)=x + x^2 + \cdots + x^{2^{m-1}}$%
- *             if %$x \in \gf{2^m}$%).
+ *             if %$x \in \gf{2^m}$%).  Since the trace is invariant under
+ *             the Frobenius automorphism (i.e., %$\Tr(x)^2 = \Tr(x)$%), it
+ *             must be an element of the base field, i.e., %$\gf{2}$%, and
+ *             we only need a single bit to represent it.
  */
 
 extern int gfreduce_trace(gfreduce */*r*/, mp */*x*/);
@@ -154,6 +157,17 @@ extern mp *gfreduce_halftrace(gfreduce */*r*/, mp */*d*/, mp */*x*/);
  *             @mp *x@ = some polynomial
  *
  * Returns:    A polynomial @y@ such that %$y^2 + y = x$%, or null.
+ *
+ * Use:                Solves quadratic equations in a field with characteristic 2.
+ *             Suppose we have an equation %$y^2 + A y + B = 0$% where
+ *             %$A \ne 0$%.  (If %$A = 0$% then %$y = \sqrt{B}$% and you
+ *             want @gfreduce_sqrt@ instead.)  Use this function to solve
+ *             %$z^2 + z = B/A^2$%; then set %$y = A z$%, since
+ *             %$y^2 + y = A^2 z^2 + A^2 z = A^2 (z^2 + z) = B$% as
+ *             required.
+ *
+ *             The two roots are %$z$% and %$z + 1$%; this function always
+ *             returns the one with zero scalar coefficient.
  */
 
 extern mp *gfreduce_quadsolve(gfreduce */*r*/, mp */*d*/, mp */*x*/);