chiark - git - mdw - rhodes/blob - rho.cc

   1 /* -*-c-*-
   2  *
   3  * Calculate discrete logs in prime-order groups
   4  *
   5  * (c) 2017 Mark Wooding
   6  */
   7
   8 /*----- Licensing notice --------------------------------------------------*
   9  *
  10  * This file is part of Rhodes, a distributed discrete-log finder.
  11  *
  12  * Rhodes is free software; you can redistribute it and/or modify
  13  * it under the terms of the GNU General Public License as published by
  14  * the Free Software Foundation; either version 2 of the License, or
  15  * (at your option) any later version.
  16  *
  17  * Rhodes 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 General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU General Public License
  23  * along with Rhodes; if not, write to the Free Software Foundation,
  24  * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  25  */
  26
  27 #include <cassert>
  28 #include <cerrno>
  29 #include <climits>
  30 #include <cstdarg>
  31 #include <cstdio>
  32 #include <cstdlib>
  33 #include <cstring>
  34
  35 #include <algorithm>
  36 #include <fstream>
  37 #include <iostream>
  38 #include <string>
  39 #include <sstream>
  40
  41 #include <sys/types.h>
  42 #include <unistd.h>
  43 #include <fcntl.h>
  44 #include <signal.h>
  45
  46 #include <NTL/GF2E.h>
  47 #include <NTL/GF2X.h>
  48 #include <NTL/GF2XFactoring.h>
  49 #include <NTL/ZZ.h>
  50 #include <NTL/ZZ_p.h>
  51
  52 static const char *prog;
  53 static int stdin_fdflags = -1;
  54
  55 static void cleanup(int sig)
  56 {
  57   if (stdin_fdflags >= 0) fcntl(0, F_SETFL, stdin_fdflags);
  58   if (sig == SIGTSTP) raise(SIGSTOP);
  59   else if (sig) { signal(sig, SIG_DFL); raise(sig); }
  60 }
  61
  62 __attribute__((noreturn))
  63 static void barf(const char *msg, int err, ...)
  64 {
  65   va_list ap;
  66
  67   va_start(ap, err);
  68   std::fprintf(stderr, "%s: ", prog);
  69   std::vfprintf(stderr, msg, ap);
  70   if (err) std::fprintf(stderr, ": %s", std::strerror(err));
  71   std::fputc('\n', stderr);
  72   va_end(ap);
  73   cleanup(0);
  74   exit(2);
  75 }
  76
  77 static void wakeup(int sig)
  78 {
  79   if (fcntl(0, F_SETFL, stdin_fdflags | O_NONBLOCK))
  80     barf("fcntl(stdin)", errno);
  81 }
  82
  83 static int intarg(const char *p, const char *what)
  84 {
  85   int e = errno;
  86   long n;
  87   char *q;
  88
  89   errno = 0;
  90   n = std::strtol(p, &q, 0);
  91   if (*q || errno || n < 0 || n > INT_MAX)
  92     barf("bad int %s `%s'", 0, what, p);
  93   errno = e;
  94   return ((int)n);
  95 }
  96
  97 static NTL::GF2X polyarg(const char *p, const char *what)
  98 {
  99   std::string s{p};
 100
 101   // Oh, this is wretched: NTL reads and writes the nibbles backwards.  Hack
 102   // hack.
 103   if (s.size() < 2 || s[0] != '0' || s[1] != 'x') goto bad;
 104   std::reverse(s.begin() + 2, s.end());
 105
 106   { NTL::GF2X a;
 107     std::istringstream ss{s};
 108     ss >> a;
 109     if (!ss) goto bad;
 110     return (a);
 111   }
 112 bad:
 113   barf("bad poly %s `%s'", 0, what, p);
 114 }
 115
 116 static unsigned long long hash(const unsigned char *p, size_t n)
 117 {
 118   size_t i;
 119   unsigned long long h = 0x6a078c523ae42f68ull;
 120
 121   for (i = 0; i < n; i++) { h += p[i]; h *= 0xbeaa913866b8d5a3ull; }
 122   return (h);
 123 }
 124
 125 static std::string showpoly(NTL::GF2X f)
 126 {
 127   std::ostringstream ss;
 128   ss << f;
 129   std::string s{ss.str()};
 130   assert(s.size() >= 2 && s[0] == '0' && s[1] == 'x');
 131   std::reverse(s.begin() + 2, s.end());
 132   return (s);
 133 }
 134
 135 static NTL::ZZ zzarg(const char *p, const char *what)
 136 {
 137   std::string s{p};
 138   std::istringstream ss{s};
 139   NTL::ZZ n;
 140   ss >> n;
 141   if (!ss) barf("bad int %s `%s'", 0, what, p);
 142   return (n);
 143 }
 144
 145 static void seed()
 146 {
 147   unsigned char b[256];
 148   FILE *fp;
 149
 150   if ((fp = std::fopen("/dev/urandom", "rb")) == 0)
 151     barf("open(/dev/urandom)", errno);
 152   if (std::fread(b, 1, sizeof(b), fp) != sizeof(b)) {
 153     barf("read(/dev/urandom)%s",
 154          ferror(fp) ? errno : 0,
 155          ferror(fp) ? "" : ": unexpected short read");
 156   }
 157   std::fclose(fp);
 158   NTL::ZZ t{NTL::ZZFromBytes(b, sizeof(b))};
 159   SetSeed(t);
 160 }
 161
 162 struct step {
 163   NTL::ZZ_p du, dv;
 164   NTL::GF2E m;
 165 };
 166
 167 struct hist {
 168   unsigned long long k;
 169   NTL::ZZ_p u, v;
 170   NTL::GF2E y;
 171 };
 172
 173 #define NSTEP 23
 174 #define NSQR 2
 175 #define NHIST 8
 176
 177 #define CHECK_NITER 65536
 178
 179 int main(int argc, char *argv[])
 180 {
 181   int i;
 182   int dpbits;
 183   unsigned char buf[4096];
 184   unsigned long long niter, dpmask;
 185   fd_set fds_in;
 186   struct timeval tv;
 187   struct sigaction sa;
 188   ssize_t n;
 189
 190   prog = argv[0];
 191   NTL::GF2X::HexOutput = 1;
 192
 193   // Collect the arguments and check that they make sense.
 194   if (argc != 7) {
 195     std::fprintf(stderr, "usage: %s DPBITS gf2x P A B L\n",
 196                  prog);
 197     std::exit(2);
 198   }
 199
 200   // Distinguished-point bits, or zero to find a full cycle.
 201   dpbits = intarg(argv[1], "dpbits");
 202   dpmask = (1ull << dpbits) - 1;
 203
 204   // Group specification.  Currently only subgroups of (GF(2)[x]/(p(x)))^*
 205   // supported.
 206   if (std::strcmp(argv[2], "gf2x") != 0)
 207     barf("unknown group kind `%s'", 0, argv[1]);
 208   NTL::GF2X p = polyarg(argv[3], "p");
 209   if (!IterIrredTest(p)) barf("not irreducible", 0);
 210   NTL::GF2E::init(p);
 211
 212   // The base and operand for the logarithm.
 213   NTL::GF2E a{to_GF2E(polyarg(argv[4], "a"))};
 214   if (a == 0 || a == 1) barf("a is trivial", 0);
 215   NTL::GF2E b{to_GF2E(polyarg(argv[5], "b"))};
 216   if (b == 0) barf("b isn't a unit", 0);
 217
 218   // The group order, which we expect to have been calculated.  This must be
 219   // prime, or we risk failing later.
 220   NTL::ZZ l{zzarg(argv[6], "l")};
 221   if (!ProbPrime(l)) barf("order isn't prime", 0);
 222   NTL::ZZ_p::init(l);
 223
 224   // Check the points at least have the same order.
 225   if (power(a, l) != 1) barf("a has wrong order", 0);
 226   if (power(b, l) != 1) barf("b has wrong order", 0);
 227
 228   // Build the table of steps.  Must do this deterministically!  The basic
 229   // structure of the walk is taken from Edlyn Teske, 1998, `Better Random
 230   // Walks for Pollard's Rho Method': about twenty random steps, together
 231   // with a couple of slots' worth of squaring.  Deviating slightly from
 232   // Teske, I'm using a prime number of slots to improve the hash
 233   // distribution.
 234   step tab[NSTEP - NSQR];
 235   SetSeed(NTL::ZZ::zero());
 236   for (i = 0; i < NSTEP - NSQR; i++) {
 237     tab[i].du = NTL::random_ZZ_p();
 238     tab[i].dv = NTL::random_ZZ_p();
 239     tab[i].m = power(a, rep(tab[i].du))*power(b, rep(tab[i].dv));
 240   }
 241
 242   // Prepare signal handlers.  We're going to make stdin be non-blocking,
 243   // which will have a persistent effect on whatever file is being used
 244   // (e.g., a terminal), so we should be careful to undo it when we exit.
 245   stdin_fdflags = fcntl(0, F_GETFL);
 246   if (stdin_fdflags < 0) barf("fcntl stdin", errno);
 247   sa.sa_handler = cleanup;
 248   sigemptyset(&sa.sa_mask);
 249   sigaddset(&sa.sa_mask, SIGTSTP);
 250   sigaddset(&sa.sa_mask, SIGCONT);
 251   sa.sa_flags = 0;
 252   if (sigaction(SIGINT, &sa, 0) ||
 253       sigaction(SIGTERM, &sa, 0) ||
 254       sigaction(SIGQUIT, &sa, 0) ||
 255       sigaction(SIGTSTP, &sa, 0))
 256     barf("sigaction", errno);
 257   sa.sa_handler = wakeup;
 258   if (sigaction(SIGCONT, &sa, 0))
 259     barf("sigaction", errno);
 260   if (fcntl(0, F_SETFL, stdin_fdflags | O_NONBLOCK))
 261     barf("fcntl stdin", errno);
 262
 263   // Now we start the random walk...
 264   seed();
 265   niter = 8ull << (dpbits ? dpbits : (NumBits(l) + 1)/2);
 266 again:
 267   NTL::ZZ_p u{NTL::random_ZZ_p()}, v{NTL::random_ZZ_p()};
 268   NTL::GF2E t = power(a, rep(u))*power(b, rep(v));
 269
 270   hist h[NHIST];
 271   unsigned o = 0;
 272   unsigned long long k = 0;
 273
 274   // Initialize the table of previous steps in the walk, if we're going to
 275   // finding a cycle.  Deviating from Teske's algorithm, initialize the table
 276   // with zeros.  If we don't do this then, in the case where b = 1, the
 277   // algorithm tends to find the cycle but can't deduce anything useful from
 278   // it.
 279   if (!dpbits)
 280     for (i = 0; i < NHIST; i++)
 281       { h[i].k = 0; h[i].y = a; h[i].u = 1; h[i].v = 0; }
 282
 283   for (;;) {
 284     if (k >= niter) goto again;
 285
 286     // Check for input on stdin.  If stdin has closed, that means that our
 287     // caller has gone away, presumably because they no longer care about
 288     // what we have to compute.  (Maybe some other job found the answer
 289     // quicker than we did.)  If stdin has data in, then read it: this will
 290     // keep a network connection alive (or notice that it's gone dead).
 291     if (!(k%CHECK_NITER)) {
 292       FD_ZERO(&fds_in); FD_SET(0, &fds_in);
 293       tv.tv_sec = 0; tv.tv_usec = 0;
 294       if (select(1, &fds_in, 0, 0, &tv) < 0)
 295         { if (errno != EINTR) barf("select", errno); }
 296       else if (FD_ISSET(0, &fds_in)) {
 297         for (;;) {
 298           n = read(0, buf, sizeof(buf));
 299           if (!n) { cleanup(0); std::exit(1); }
 300           else if (n > 0) continue;
 301           else if (errno == EAGAIN) break;
 302           else if (errno != EINTR) barf("read stdin", errno);
 303         }
 304       }
 305     }
 306     k++;
 307
 308     size_t nb = NumBytes(rep(t));
 309     BytesFromGF2X(buf, rep(t), nb);
 310     unsigned long long hh = hash(buf, nb);
 311
 312     if (dpbits) {
 313       // If we're walking to a distinguished point (the algorithm of Paul van
 314       // Oorschot and Michael Wiener, 1994, `Parallel Collision Search with
 315       // Application to Hash Functions and Discrete Logarithms') then check
 316       // the hash to see if enough low-order bits are zero.
 317       if (!(hh&dpmask)) {
 318         std::cout << u << " " << v << " " << showpoly(rep(t)) << " "
 319                   << k << std::endl;
 320         goto done;
 321       }
 322     } else {
 323       // Otherwise, check for a cycle.  We use the algorithm is from Edlyn
 324       // Teske, 1998, `A Space Efficient Algorithm for Group Structure
 325       // Computation', rather than the usual one by Floyd, which requires
 326       // about three times as much computation.
 327
 328       // Check for a match.
 329       for (i = 0; i < NHIST; i++) {
 330         if (t == h[i].y) {
 331           if (h[i].u == u || h[i].v == v) goto again;
 332           std::cout << (h[i].u - u)/(v - h[i].v) << " " << k << std::endl;
 333           goto done;
 334         }
 335       }
 336
 337       // Update the table of earlier points if necessary.
 338       if (k > 3*h[o].k) {
 339         h[o].u = u; h[o].v = v; h[o].y = t; h[o].k = k;
 340         o = (o + 1)%NHIST;
 341       }
 342     }
 343
 344     // Take a step in the random walk.
 345     i = hh%NSTEP;
 346     if (i >= NSTEP - NSQR) { mul(u, u, 2); mul(v, v, 2); sqr(t, t);; }
 347     else { add(u, u, tab[i].du); add(v, v, tab[i].dv); mul(t, t, tab[i].m); }
 348   }
 349
 350 done:
 351   cleanup(0);
 352   return (0);
 353 }