ec-field-test.c: Make the field-element type use internal format.
[secnet] / fgoldi.c
index 6a8d35b..08e6592 100644 (file)
--- a/fgoldi.c
+++ b/fgoldi.c
@@ -1,41 +1,3 @@
-/*
- * fgoldi.c: arithmetic modulo 2^448 - 2^224 - 1
- */
-/*
- * This file is Free Software.  It has been modified to as part of its
- * incorporation into secnet.
- *
- * Copyright 2017 Mark Wooding
- *
- * You may redistribute this file and/or modify it under the terms of
- * the permissive licence shown below.
- *
- * You may redistribute secnet as a whole and/or modify it under the
- * terms of the GNU General Public License as published by the Free
- * Software Foundation; either version 3, or (at your option) any
- * later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, see
- * https://www.gnu.org/licenses/gpl.html.
- */
-/*
- * Imported from Catacomb, and modified for Secnet (2017-04-30):
- *
- *   * Use `fake-mLib-bits.h' in place of the real <mLib/bits.h>.
- *
- *   * Remove the 16/32-bit implementation, since C99 always has 64-bit
- *     arithmetic.
- *
- *   * Remove the test rig code: a replacement is in a separate source file.
- *
- * The file's original comment headers are preserved below.
- */
 /* -*-c-*-
  *
  * Arithmetic in the Goldilocks field GF(2^448 - 2^224 - 1)
@@ -45,7 +7,26 @@
 
 /*----- Licensing notice --------------------------------------------------*
  *
- * This file is part of Catacomb.
+ * This file is part of secnet.
+ * See README for full list of copyright holders.
+ *
+ * secnet is free software; you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version d of the License, or
+ * (at your option) any later version.
+ *
+ * secnet is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * version 3 along with secnet; if not, see
+ * https://www.gnu.org/licenses/gpl.html.
+ *
+ * This file was originally part of Catacomb, but has been automatically
+ * modified for incorporation into secnet: see `import-catacomb-crypto'
+ * for details.
  *
  * Catacomb is free software; you can redistribute it and/or modify
  * it under the terms of the GNU Library General Public License as
  * (hence the name).
  */
 
+typedef fgoldi_piece piece;
+
 /* We represent an element of GF(p) as 16 28-bit signed integer pieces x_i:
  * x = SUM_{0<=i<16} x_i 2^(28i).
  */
 
-typedef  int32  piece; typedef  int64  dblpiece;
+                      typedef  int64  dblpiece;
 typedef uint32 upiece; typedef uint64 udblpiece;
 #define PIECEWD(i) 28
 #define NPIECE 16
 #define P p28
 
-#define B28 0x10000000u
 #define B27 0x08000000u
 #define M28 0x0fffffffu
-#define M27 0x07ffffffu
 #define M32 0xffffffffu
 
 /*----- Debugging machinery -----------------------------------------------*/
@@ -135,6 +116,7 @@ DEF_FDUMP(fdump, piece, PIECEWD, NPIECE, 56, get_pgoldi())
 
 void fgoldi_load(fgoldi *z, const octet xv[56])
 {
+
   unsigned i;
   uint32 xw[14];
   piece b, c;
@@ -183,6 +165,7 @@ void fgoldi_load(fgoldi *z, const octet xv[56])
 
 void fgoldi_store(octet zv[56], const fgoldi *x)
 {
+
   piece y[NPIECE], yy[NPIECE], c, d;
   uint32 u, v;
   mask32 m;
@@ -302,8 +285,72 @@ void fgoldi_sub(fgoldi *z, const fgoldi *x, const fgoldi *y)
   for (i = 0; i < NPIECE; i++) z->P[i] = x->P[i] - y->P[i];
 }
 
+/* --- @fgoldi_neg@ --- *
+ *
+ * Arguments:  @fgoldi *z@ = where to put the result (may alias @x@)
+ *             @const fgoldi *x@ = an operand
+ *
+ * Returns:    ---
+ *
+ * Use:                Set @z = -x@.
+ */
+
+void fgoldi_neg(fgoldi *z, const fgoldi *x)
+{
+  unsigned i;
+  for (i = 0; i < NPIECE; i++) z->P[i] = -x->P[i];
+}
+
 /*----- Constant-time utilities -------------------------------------------*/
 
+/* --- @fgoldi_pick2@ --- *
+ *
+ * Arguments:  @fgoldi *z@ = where to put the result (may alias @x@ or @y@)
+ *             @const fgoldi *x, *y@ = two operands
+ *             @uint32 m@ = a mask
+ *
+ * Returns:    ---
+ *
+ * Use:                If @m@ is zero, set @z = y@; if @m@ is all-bits-set, then set
+ *             @z = x@.  If @m@ has some other value, then scramble @z@ in
+ *             an unhelpful way.
+ */
+
+void fgoldi_pick2(fgoldi *z, const fgoldi *x, const fgoldi *y, uint32 m)
+{
+  mask32 mm = FIX_MASK32(m);
+  unsigned i;
+  for (i = 0; i < NPIECE; i++) z->P[i] = PICK2(x->P[i], y->P[i], mm);
+}
+
+/* --- @fgoldi_pickn@ --- *
+ *
+ * Arguments:  @fgoldi *z@ = where to put the result
+ *             @const fgoldi *v@ = a table of entries
+ *             @size_t n@ = the number of entries in @v@
+ *             @size_t i@ = an index
+ *
+ * Returns:    ---
+ *
+ * Use:                If @0 <= i < n < 32@ then set @z = v[i]@.  If @n >= 32@ then
+ *             do something unhelpful; otherwise, if @i >= n@ then set @z@
+ *             to zero.
+ */
+
+void fgoldi_pickn(fgoldi *z, const fgoldi *v, size_t n, size_t i)
+{
+  uint32 b = (uint32)1 << (31 - i);
+  mask32 m;
+  unsigned j;
+
+  for (j = 0; j < NPIECE; j++) z->P[j] = 0;
+  while (n--) {
+    m = SIGN(b);
+    for (j = 0; j < NPIECE; j++) CONDPICK(z->P[j], v->P[j], m);
+    v++; b <<= 1;
+  }
+}
+
 /* --- @fgoldi_condswap@ --- *
  *
  * Arguments:  @fgoldi *x, *y@ = two operands
@@ -324,6 +371,31 @@ void fgoldi_condswap(fgoldi *x, fgoldi *y, uint32 m)
   for (i = 0; i < NPIECE; i++) CONDSWAP(x->P[i], y->P[i], mm);
 }
 
+/* --- @fgoldi_condneg@ --- *
+ *
+ * Arguments:  @fgoldi *z@ = where to put the result (may alias @x@)
+ *             @const fgoldi *x@ = an operand
+ *             @uint32 m@ = a mask
+ *
+ * Returns:    ---
+ *
+ * Use:                If @m@ is zero, set @z = x@; if @m@ is all-bits-set, then set
+ *             @z = -x@.  If @m@ has some other value then scramble @z@ in
+ *             an unhelpful way.
+ */
+
+void fgoldi_condneg(fgoldi *z, const fgoldi *x, uint32 m)
+{
+  mask32 m_xor = FIX_MASK32(m);
+  piece m_add = m&1;
+# define CONDNEG(x) (((x) ^ m_xor) + m_add)
+
+  unsigned i;
+  for (i = 0; i < NPIECE; i++) z->P[i] = CONDNEG(x->P[i]);
+
+#undef CONDNEG
+}
+
 /*----- Multiplication ----------------------------------------------------*/
 
 /* Let B = 2^63 - 1 be the largest value such that +B and -B can be
@@ -473,6 +545,7 @@ void fgoldi_mul(fgoldi *z, const fgoldi *x, const fgoldi *y)
 
 void fgoldi_sqr(fgoldi *z, const fgoldi *x)
 {
+
   dblpiece zz[NPIECE], u[NPIECE];
   piece ab[NPIECE];
   const piece *a = x->P + NPIECE/2, *b = x->P;
@@ -543,7 +616,6 @@ void fgoldi_sqr(fgoldi *z, const fgoldi *x)
   /* Finally, carrying. */
   for (i = 0; i < 2; i++) CARRY_REDUCE(zz, zz);
   for (i = 0; i < NPIECE; i++) z->P[i] = zz[i];
-
 }
 
 /*----- More advanced operations ------------------------------------------*/
@@ -603,4 +675,86 @@ void fgoldi_inv(fgoldi *z, const fgoldi *x)
 #undef SQRN
 }
 
+/* --- @fgoldi_quosqrt@ --- *
+ *
+ * Arguments:  @fgoldi *z@ = where to put the result (may alias @x@ or @y@)
+ *             @const fgoldi *x, *y@ = two operands
+ *
+ * Returns:    Zero if successful, @-1@ if %$x/y$% is not a square.
+ *
+ * Use:                Stores in @z@ the one of the square roots %$\pm\sqrt{x/y}$%.
+ *             If %$x = y = 0% then the result is zero; if %$y = 0$% but %$x
+ *             \ne 0$% then the operation fails.  If you wanted a specific
+ *             square root then you'll have to pick it yourself.
+ */
+
+int fgoldi_quosqrt(fgoldi *z, const fgoldi *x, const fgoldi *y)
+{
+  fgoldi t, u, v;
+  octet xb[56], b0[56];
+  int32 rc = -1;
+  mask32 m;
+  unsigned i;
+
+#define SQRN(z, x, n) do {                                             \
+  fgoldi_sqr((z), (x));                                                        \
+  for (i = 1; i < (n); i++) fgoldi_sqr((z), (z));                      \
+} while (0)
+
+  /* This is, fortunately, significantly easier than the equivalent problem
+   * in GF(2^255 - 19), since p == 3 (mod 4).
+   *
+   * If x/y is square, then so is v = y^2 x/y = x y, and therefore u has
+   * order r = (p - 1)/2.  Let w = v^{(p-3)/4}.  Then w^2 = v^{(p-3)/2} =
+   * u^{r-1} = 1/v = 1/x y.  Clearly, then, (x w)^2 = x^2/x y = x/y, so x w
+   * is one of the square roots we seek.
+   *
+   * The addition chain, then, is a prefix of the previous one.
+   */
+  fgoldi_mul(&v, x, y);
+
+  fgoldi_sqr(&u, &v);                  /*    1 | 2 */
+  fgoldi_mul(&t, &u, &v);              /*    2 | 3 */
+  SQRN(&u, &t, 2);                     /*    4 | 12 */
+  fgoldi_mul(&t, &u, &t);              /*    5 | 15 */
+  SQRN(&u, &t, 4);                     /*    9 | 240 */
+  fgoldi_mul(&u, &u, &t);              /*   10 | 255 = 2^8 - 1 */
+  SQRN(&u, &u, 4);                     /*   14 | 2^12 - 16 */
+  fgoldi_mul(&t, &u, &t);              /*   15 | 2^12 - 1 */
+  SQRN(&u, &t, 12);                    /*   27 | 2^24 - 2^12 */
+  fgoldi_mul(&u, &u, &t);              /*   28 | 2^24 - 1 */
+  SQRN(&u, &u, 12);                    /*   40 | 2^36 - 2^12 */
+  fgoldi_mul(&t, &u, &t);              /*   41 | 2^36 - 1 */
+  fgoldi_sqr(&t, &t);                  /*   42 | 2^37 - 2 */
+  fgoldi_mul(&t, &t, &v);              /*   43 | 2^37 - 1 */
+  SQRN(&u, &t, 37);                    /*   80 | 2^74 - 2^37 */
+  fgoldi_mul(&u, &u, &t);              /*   81 | 2^74 - 1 */
+  SQRN(&u, &u, 37);                    /*  118 | 2^111 - 2^37 */
+  fgoldi_mul(&t, &u, &t);              /*  119 | 2^111 - 1 */
+  SQRN(&u, &t, 111);                   /*  230 | 2^222 - 2^111 */
+  fgoldi_mul(&t, &u, &t);              /*  231 | 2^222 - 1 */
+  fgoldi_sqr(&u, &t);                  /*  232 | 2^223 - 2 */
+  fgoldi_mul(&u, &u, &v);              /*  233 | 2^223 - 1 */
+  SQRN(&u, &u, 223);                   /*  456 | 2^446 - 2^223 */
+  fgoldi_mul(&t, &u, &t);              /*  457 | 2^446 - 2^222 - 1 */
+
+#undef SQRN
+
+  /* Now we must decide whether the answer was right.  We should have z^2 =
+   * x/y, so y z^2 = x.
+   *
+   * The easiest way to compare is to encode.  This isn't as wasteful as it
+   * sounds: the hard part is normalizing the representations, which we have
+   * to do anyway.
+   */
+  fgoldi_mul(z, x, &t);
+  fgoldi_sqr(&t, z);
+  fgoldi_mul(&t, &t, y);
+  fgoldi_store(xb, x);
+  fgoldi_store(b0, &t);
+  m = -consttime_memeq(xb, b0, 56);
+  rc = PICK2(0, rc, m);
+  return (rc);
+}
+
 /*----- That's all, folks -------------------------------------------------*/