Add commentary and licence notices.

[rhodes] / rho.cc

/* -*-c-*-
 *
 * Calculate discrete logs in prime-order groups
 *
 * (c) 2017 Mark Wooding
 */

/*----- Licensing notice --------------------------------------------------*
 *
 * This file is part of Rhodes, a distributed discrete-log finder.
 *
 * Rhodes 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 2 of the License, or
 * (at your option) any later version.
 *
 * Rhodes 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 Rhodes; if not, write to the Free Software Foundation,
 * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 */

#include <cassert>
#include <cerrno>
#include <climits>
#include <cstdarg>
#include <cstdio>
#include <cstdlib>
#include <cstring>

#include <algorithm>
#include <fstream>
#include <iostream>
#include <string>
#include <sstream>

#include <sys/types.h>
#include <unistd.h>
#include <fcntl.h>
#include <signal.h>

#include <NTL/GF2E.h>
#include <NTL/GF2X.h>
#include <NTL/GF2XFactoring.h>
#include <NTL/ZZ.h>
#include <NTL/ZZ_p.h>

static const char *prog;
static int stdin_fdflags = -1;

static void cleanup(int sig)
{
  if (stdin_fdflags >= 0) fcntl(0, F_SETFL, stdin_fdflags);
  if (sig == SIGTSTP) raise(SIGSTOP);
  else if (sig) { signal(sig, SIG_DFL); raise(sig); }
}

__attribute__((noreturn))
static void barf(const char *msg, int err, ...)
{
  va_list ap;

  va_start(ap, err);
  std::fprintf(stderr, "%s: ", prog);
  std::vfprintf(stderr, msg, ap);
  if (err) std::fprintf(stderr, ": %s", std::strerror(err));
  std::fputc('\n', stderr);
  va_end(ap);
  cleanup(0);
  exit(2);
}

static void wakeup(int sig)
{
  if (fcntl(0, F_SETFL, stdin_fdflags | O_NONBLOCK))
    barf("fcntl(stdin)", errno);
}

static int intarg(const char *p, const char *what)
{
  int e = errno;
  long n;
  char *q;

  errno = 0;
  n = std::strtol(p, &q, 0);
  if (*q || errno || n < 0 || n > INT_MAX)
    barf("bad int %s `%s'", 0, what, p);
  errno = e;
  return ((int)n);
}

static NTL::GF2X polyarg(const char *p, const char *what)
{
  std::string s{p};

  // Oh, this is wretched: NTL reads and writes the nibbles backwards.  Hack
  // hack.
  if (s.size() < 2 || s[0] != '0' || s[1] != 'x') goto bad;
  std::reverse(s.begin() + 2, s.end());

  { NTL::GF2X a;
    std::istringstream ss{s};
    ss >> a;
    if (!ss) goto bad;
    return (a);
  }
bad:
  barf("bad poly %s `%s'", 0, what, p);
}

static unsigned long long hash(const unsigned char *p, size_t n)
{
  size_t i;
  unsigned long long h = 0x6a078c523ae42f68ull;

  for (i = 0; i < n; i++) { h += p[i]; h *= 0xbeaa913866b8d5a3ull; }
  return (h);
}

static std::string showpoly(NTL::GF2X f)
{
  std::ostringstream ss;
  ss << f;
  std::string s{ss.str()};
  assert(s.size() >= 2 && s[0] == '0' && s[1] == 'x');
  std::reverse(s.begin() + 2, s.end());
  return (s);
}

static NTL::ZZ zzarg(const char *p, const char *what)
{
  std::string s{p};
  std::istringstream ss{s};
  NTL::ZZ n;
  ss >> n;
  if (!ss) barf("bad int %s `%s'", 0, what, p);
  return (n);
}

static void seed()
{
  unsigned char b[256];
  FILE *fp;

  if ((fp = std::fopen("/dev/urandom", "rb")) == 0)
    barf("open(/dev/urandom)", errno);
  if (std::fread(b, 1, sizeof(b), fp) != sizeof(b)) {
    barf("read(/dev/urandom)%s",
	 ferror(fp) ? errno : 0,
	 ferror(fp) ? "" : ": unexpected short read");
  }
  std::fclose(fp);
  NTL::ZZ t{NTL::ZZFromBytes(b, sizeof(b))};
  SetSeed(t);
}

struct step {
  NTL::ZZ_p du, dv;
  NTL::GF2E m;
};

struct hist {
  unsigned long long k;
  NTL::ZZ_p u, v;
  NTL::GF2E y;
};

#define NSTEP 23
#define NSQR 2
#define NHIST 8

#define CHECK_NITER 65536

int main(int argc, char *argv[])
{
  int i;
  int dpbits;
  unsigned char buf[4096];
  unsigned long long niter, dpmask;
  fd_set fds_in;
  struct timeval tv;
  struct sigaction sa;
  ssize_t n;

  prog = argv[0];
  NTL::GF2X::HexOutput = 1;

  // Collect the arguments and check that they make sense.
  if (argc != 7) {
    std::fprintf(stderr, "usage: %s DPBITS gf2x P A B L\n",
		 prog);
    std::exit(2);
  }

  // Distinguished-point bits, or zero to find a full cycle.
  dpbits = intarg(argv[1], "dpbits");
  dpmask = (1ull << dpbits) - 1;

  // Group specification.  Currently only subgroups of (GF(2)[x]/(p(x)))^*
  // supported.
  if (std::strcmp(argv[2], "gf2x") != 0)
    barf("unknown group kind `%s'", 0, argv[1]);
  NTL::GF2X p = polyarg(argv[3], "p");
  if (!IterIrredTest(p)) barf("not irreducible", 0);
  NTL::GF2E::init(p);

  // The base and operand for the logarithm.
  NTL::GF2E a{to_GF2E(polyarg(argv[4], "a"))};
  if (a == 0 || a == 1) barf("a is trivial", 0);
  NTL::GF2E b{to_GF2E(polyarg(argv[5], "b"))};
  if (b == 0) barf("b isn't a unit", 0);

  // The group order, which we expect to have been calculated.  This must be
  // prime, or we risk failing later.
  NTL::ZZ l{zzarg(argv[6], "l")};
  if (!ProbPrime(l)) barf("order isn't prime", 0);
  NTL::ZZ_p::init(l);

  // Check the points at least have the same order.
  if (power(a, l) != 1) barf("a has wrong order", 0);
  if (power(b, l) != 1) barf("b has wrong order", 0);

  // Build the table of steps.  Must do this deterministically!  The basic
  // structure of the walk is taken from Edlyn Teske, 1998, `Better Random
  // Walks for Pollard's Rho Method': about twenty random steps, together
  // with a couple of slots' worth of squaring.  Deviating slightly from
  // Teske, I'm using a prime number of slots to improve the hash
  // distribution.
  step tab[NSTEP - NSQR];
  SetSeed(NTL::ZZ::zero());
  for (i = 0; i < NSTEP - NSQR; i++) {
    tab[i].du = NTL::random_ZZ_p();
    tab[i].dv = NTL::random_ZZ_p();
    tab[i].m = power(a, rep(tab[i].du))*power(b, rep(tab[i].dv));
  }

  // Prepare signal handlers.  We're going to make stdin be non-blocking,
  // which will have a persistent effect on whatever file is being used
  // (e.g., a terminal), so we should be careful to undo it when we exit.
  stdin_fdflags = fcntl(0, F_GETFL);
  if (stdin_fdflags < 0) barf("fcntl stdin", errno);
  sa.sa_handler = cleanup;
  sigemptyset(&sa.sa_mask);
  sigaddset(&sa.sa_mask, SIGTSTP);
  sigaddset(&sa.sa_mask, SIGCONT);
  sa.sa_flags = 0;
  if (sigaction(SIGINT, &sa, 0) ||
      sigaction(SIGTERM, &sa, 0) ||
      sigaction(SIGQUIT, &sa, 0) ||
      sigaction(SIGTSTP, &sa, 0))
    barf("sigaction", errno);
  sa.sa_handler = wakeup;
  if (sigaction(SIGCONT, &sa, 0))
    barf("sigaction", errno);
  if (fcntl(0, F_SETFL, stdin_fdflags | O_NONBLOCK))
    barf("fcntl stdin", errno);

  // Now we start the random walk...
  seed();
  niter = 8ull << (dpbits ? dpbits : (NumBits(l) + 1)/2);
again:
  NTL::ZZ_p u{NTL::random_ZZ_p()}, v{NTL::random_ZZ_p()};
  NTL::GF2E t = power(a, rep(u))*power(b, rep(v));

  hist h[NHIST];
  unsigned o = 0;
  unsigned long long k = 0;

  // Initialize the table of previous steps in the walk, if we're going to
  // finding a cycle.  Deviating from Teske's algorithm, initialize the table
  // with zeros.  If we don't do this then, in the case where b = 1, the
  // algorithm tends to find the cycle but can't deduce anything useful from
  // it.
  if (!dpbits)
    for (i = 0; i < NHIST; i++)
      { h[i].k = 0; h[i].y = a; h[i].u = 1; h[i].v = 0; }

  for (;;) {
    if (k >= niter) goto again;

    // Check for input on stdin.  If stdin has closed, that means that our
    // caller has gone away, presumably because they no longer care about
    // what we have to compute.  (Maybe some other job found the answer
    // quicker than we did.)  If stdin has data in, then read it: this will
    // keep a network connection alive (or notice that it's gone dead).
    if (!(k%CHECK_NITER)) {
      FD_ZERO(&fds_in); FD_SET(0, &fds_in);
      tv.tv_sec = 0; tv.tv_usec = 0;
      if (select(1, &fds_in, 0, 0, &tv) < 0)
	{ if (errno != EINTR) barf("select", errno); }
      else if (FD_ISSET(0, &fds_in)) {
	for (;;) {
	  n = read(0, buf, sizeof(buf));
	  if (!n) { cleanup(0); std::exit(1); }
	  else if (n > 0) continue;
	  else if (errno == EAGAIN) break;
	  else if (errno != EINTR) barf("read stdin", errno);
	}
      }
    }
    k++;

    size_t nb = NumBytes(rep(t));
    BytesFromGF2X(buf, rep(t), nb);
    unsigned long long hh = hash(buf, nb);

    if (dpbits) {
      // If we're walking to a distinguished point (the algorithm of Paul van
      // Oorschot and Michael Wiener, 1994, `Parallel Collision Search with
      // Application to Hash Functions and Discrete Logarithms') then check
      // the hash to see if enough low-order bits are zero.
      if (!(hh&dpmask)) {
	std::cout << u << " " << v << " " << showpoly(rep(t)) << " "
		  << k << std::endl;
	goto done;
      }
    } else {
      // Otherwise, check for a cycle.  We use the algorithm is from Edlyn
      // Teske, 1998, `A Space Efficient Algorithm for Group Structure
      // Computation', rather than the usual one by Floyd, which requires
      // about three times as much computation.

      // Check for a match.
      for (i = 0; i < NHIST; i++) {
	if (t == h[i].y) {
	  if (h[i].u == u || h[i].v == v) goto again;
	  std::cout << (h[i].u - u)/(v - h[i].v) << " " << k << std::endl;
	  goto done;
	}
      }

      // Update the table of earlier points if necessary.
      if (k > 3*h[o].k) {
	h[o].u = u; h[o].v = v; h[o].y = t; h[o].k = k;
	o = (o + 1)%NHIST;
      }
    }

    // Take a step in the random walk.
    i = hh%NSTEP;
    if (i >= NSTEP - NSQR) { mul(u, u, 2); mul(v, v, 2); sqr(t, t);; }
    else { add(u, u, tab[i].du); add(v, v, tab[i].dv); mul(t, t, tab[i].m); }
  }

done:
  cleanup(0);
  return (0);
}