chiark - git - mdw - catacomb/blob - math/qfarith.h

   1 /* -*-c-*-
   2  *
   3  * Utilities for quick field arithmetic
   4  *
   5  * (c) 2017 Straylight/Edgeware
   6  */
   7
   8 /*----- Licensing notice --------------------------------------------------*
   9  *
  10  * This file is part of Catacomb.
  11  *
  12  * Catacomb is free software; you can redistribute it and/or modify
  13  * it under the terms of the GNU Library General Public License as
  14  * published by the Free Software Foundation; either version 2 of the
  15  * License, or (at your option) any later version.
  16  *
  17  * Catacomb is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  20  * GNU Library General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Library General Public
  23  * License along with Catacomb; if not, write to the Free
  24  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  25  * MA 02111-1307, USA.
  26  */
  27
  28 #ifndef CATACOMB_QFARITH_H
  29 #define CATACOMB_QFARITH_H
  30
  31 #ifdef __cplusplus
  32   extern "C" {
  33 #endif
  34
  35 /*----- Header files ------------------------------------------------------*/
  36
  37 #include <limits.h>
  38
  39 #include <mLib/bits.h>
  40
  41 /*----- Signed integer types ----------------------------------------------*/
  42
  43 /* See if we can find a suitable 64-bit or wider type.  Don't bother if we
  44  * don't have a corresponding unsigned type, because we really need both.
  45  */
  46 #ifdef HAVE_UINT64
  47 #  if INT_MAX >> 31 == 0xffffffff
  48 #    define HAVE_INT64
  49      typedef int int64;
  50 #  elif LONG_MAX >> 31 == 0xffffffff
  51 #    define HAVE_INT64
  52      typedef long int64;
  53 #  elif defined(LLONG_MAX)
  54 #    define HAVE_INT64
  55      MLIB_BITS_EXTENSION typedef long long int64;
  56 #  endif
  57 #endif
  58
  59 /* Choose suitable 32- and 16-bit types. */
  60 #if INT_MAX >= 0x7fffffff
  61    typedef int int32;
  62 #else
  63    typedef long int32;
  64 #endif
  65
  66 typedef short int16;
  67
  68 /*----- General bit-hacking utilities -------------------------------------*/
  69
  70 /* Individual bits, and masks for low bits. */
  71 #define BIT(n) (1ul << (n))
  72 #define MASK(n) (BIT(n) - 1)
  73
  74 /* Arithmetic right shift.  If X is a value of type TY, and N is a
  75  * nonnegative integer, then return the value of X shifted right by N bits;
  76  * alternatively, this is floor(X/2^N).
  77  *
  78  * GCC manages to compile this into a simple shift operation if one is
  79  * available, but it's correct for all C targets.
  80  */
  81 #define ASR(ty, x, n) (((x) - (ty)((u##ty)(x)&MASK(n)))/(ty)BIT(n))
  82
  83 /*----- Constant-time utilities -------------------------------------------*/
  84
  85 /* The following have better implementations on a two's complement target. */
  86 #ifdef NEG_TWOC
  87
  88   /* If we have two's complement arithmetic then masks are signed; this
  89    * avoids a switch to unsigned representation, with the consequent problem
  90    * of overflow when we convert back.
  91    */
  92   typedef int32 mask32;
  93
  94   /* Convert an unsigned mask M into a `mask32'.  This is a hairy-looking
  95    * no-op on many targets, but, given that we have two's complement
  96    * integers, it is free of arithmetic overflow.
  97    */
  98 # define FIX_MASK32(m) \
  99     ((mask32)((m)&0x7fffffffu) + (-(mask32)0x7fffffff - 1)*(((m) >> 31)&1u))
 100
 101   /* If Z is zero and M has its low 32 bits set, then copy (at least) the low
 102    * 32 bits of X to Z; if M is zero, do nothing.  Otherwise, scramble Z
 103    * unhelpfully.
 104    */
 105 # define CONDPICK(z, x, m) do { (z) |= (x)&(m); } while (0)
 106
 107   /* If M has its low 32 bits set, then return (at least) the low 32 bits of
 108    * X in Z; if M is zero, then return (at least) the low 32 bits of Y in Z.
 109    * Otherwise, return an unhelful result.
 110    */
 111 # define PICK2(x, y, m) (((x)&(m)) | ((y)&~(m)))
 112
 113   /* If M has its low 32 bits set then swap (at least) the low 32 bits of X
 114    * and Y; if M is zero, then do nothing.  Otherwise, scramble X and Y
 115    * unhelpfully.
 116    */
 117 # define CONDSWAP(x, y, m)                                              \
 118     do { mask32 t_ = ((x) ^ (y))&(m); (x) ^= t_; (y) ^= t_; } while (0)
 119 #else
 120
 121   /* We don't have two's complement arithmetic.  We can't use bithacking at
 122    * all: if we try to hack on the bits of signed numbers we'll come unstuck
 123    * when we hit the other representation of zero; and if we switch to
 124    * unsigned arithmetic then we'll have overflow when we try to convert a
 125    * negative number back.  So fall back to simple arithmetic.
 126    */
 127   typedef uint32 mask32;
 128
 129   /* Convert an unsigned mask M into a `mask32'.  Our masks are unsigned, so
 130    * this does nothing.
 131    */
 132 # define FIX_MASK32(m) ((mask32)(m))
 133
 134   /* If Z is zero and M has its low 32 bits set, then copy (at least) the low
 135    * 32 bits of X to Z; if M is zero, do nothing.  Otherwise, scramble Z
 136    * unhelpfully.
 137    */
 138 # define CONDPICK(z, x, m)                                              \
 139     do { (z) += (x)*(int)((unsigned)(m)&1u); } while (0)
 140
 141   /* If M has its low 32 bits set, then return (at least) the low 32 bits of
 142    * X in Z; if M is zero, then return (at least) the low 32 bits of Y in Z.
 143    * Otherwise, return an unhelful result.
 144    */
 145 # define PICK2(x, y, m)                                                 \
 146     ((x)*(int)((unsigned)(m)&1u) + (y)*(int)(1 - ((unsigned)(m)&1u)))
 147
 148   /* If M has its low 32 bits set then swap (at least) the low 32 bits of X
 149    * and Y; if M is zero, then do nothing.  Otherwise, scramble X and Y
 150    * unhelpfully.
 151    */
 152 # define CONDSWAP(x, y, m) do {                                         \
 153     int32 x_ = PICK2((y), (x), (m)), y_ = PICK2((x), (y), (m));         \
 154     (x) = x_; (y) = y_;                                                 \
 155   } while (0)
 156 #endif
 157
 158 /* Return zero if bit 31 of X is clear, or a mask with (at least) the low 32
 159  * bits set if bit 31 of X is set.
 160  */
 161 #define SIGN(x) (-(mask32)(((uint32)(x) >> 31)&1))
 162
 163 /* Return zero if X is zero, or a mask with (at least) the low 32 bits set if
 164  * X is nonzero.
 165  */
 166 #define NONZEROP(x) SIGN((U32(x) >> 1) - U32(x))
 167
 168 /*----- Debugging utilities -----------------------------------------------*/
 169
 170 /* Define a debugging function DUMPFN, which will dump an integer represented
 171  * modulo M.  The integer is represented as a vector of NPIECE pieces of type
 172  * PIECETY.  The pieces are assembled at possibly irregular offsets: piece i
 173  * logically has width PIECEWD(i), but may overhang the next piece.  The
 174  * pieces may be signed.  GETMOD is an expression which calculates and
 175  * returns the value of M, as an `mp *'.
 176  *
 177  * The generated function writes the value of such an integer X to the stream
 178  * FP, labelled with the string WHAT.
 179  *
 180  * The definition assumes that <stdio.h>, <catacomb/mp.h>, and
 181  * <catacomb/mptext.h> have been included.
 182  */
 183 #define DEF_FDUMP(dumpfn, piecety, piecewd, npiece, noctet, getmod)     \
 184   static void dumpfn(FILE *fp, const char *what, const piecety *x)      \
 185   {                                                                     \
 186     mpw w;                                                              \
 187     mp m, *y = MP_ZERO, *t = MP_NEW, *p;                                \
 188     octet b[noctet];                                                    \
 189     unsigned i, o;                                                      \
 190                                                                         \
 191     p = getmod;                                                         \
 192     mp_build(&m, &w, &w + 1);                                           \
 193     for (i = o = 0; i < npiece; i++) {                                  \
 194       if (x[i] >= 0) { w = x[i]; m.f &= ~MP_NEG; }                      \
 195       else { w = -x[i]; m.f |= MP_NEG; }                                \
 196       t = mp_lsl(t, &m, o);                                             \
 197       y = mp_add(y, y, t);                                              \
 198       o += piecewd(i);                                                  \
 199     }                                                                   \
 200                                                                         \
 201     fprintf(fp, "%s = <", what);                                        \
 202     for (i = 0; i < npiece; i++) {                                      \
 203       if (i) fputs(", ", fp);                                           \
 204       fprintf(fp, "%ld", (long)x[i]);                                   \
 205     }                                                                   \
 206     fputs(">\n\t= ", fp);                                               \
 207     mp_writefile(y, fp, 10);                                            \
 208     fputs("\n\t== ", fp);                                               \
 209     mp_div(0, &y, y, p);                                                \
 210     mp_writefile(y, fp, 10);                                            \
 211     fputs("\n\t= 0x", fp);                                              \
 212     mp_writefile(y, fp, 16);                                            \
 213     fputs(" (mod 2^255 - 19)\n\t= [", fp);                              \
 214     mp_storel(y, b, sizeof(b));                                         \
 215     for (i = 0; i < 32; i++) {                                          \
 216       if (i && !(i%4)) fputc(':', fp);                                  \
 217       fprintf(fp, "%02x", b[i]);                                        \
 218     }                                                                   \
 219     fputs("]\n", fp);                                                   \
 220     mp_drop(y); mp_drop(p); mp_drop(t);                                 \
 221   }
 222
 223 /*----- That's all, folks -------------------------------------------------*/
 224
 225 #ifdef __cplusplus
 226   }
 227 #endif
 228
 229 #endif