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

   1 #include <cassert>
   2 #include <cerrno>
   3 #include <climits>
   4 #include <cstdarg>
   5 #include <cstdio>
   6 #include <cstdlib>
   7 #include <cstring>
   8
   9 #include <algorithm>
  10 #include <fstream>
  11 #include <iostream>
  12 #include <string>
  13 #include <sstream>
  14
  15 #include <sys/types.h>
  16 #include <unistd.h>
  17 #include <fcntl.h>
  18 #include <signal.h>
  19
  20 #include <NTL/GF2E.h>
  21 #include <NTL/GF2X.h>
  22 #include <NTL/GF2XFactoring.h>
  23 #include <NTL/ZZ.h>
  24 #include <NTL/ZZ_p.h>
  25
  26 static const char *prog;
  27 static int stdin_fdflags = -1;
  28
  29 static void cleanup(int sig)
  30 {
  31   if (stdin_fdflags >= 0) fcntl(0, F_SETFL, stdin_fdflags);
  32   if (sig == SIGTSTP) raise(SIGSTOP);
  33   else if (sig) { signal(sig, SIG_DFL); raise(sig); }
  34 }
  35
  36 __attribute__((noreturn))
  37 static void barf(const char *msg, int err, ...)
  38 {
  39   va_list ap;
  40
  41   va_start(ap, err);
  42   std::fprintf(stderr, "%s: ", prog);
  43   std::vfprintf(stderr, msg, ap);
  44   if (err) std::fprintf(stderr, ": %s", std::strerror(err));
  45   std::fputc('\n', stderr);
  46   va_end(ap);
  47   cleanup(0);
  48   exit(2);
  49 }
  50
  51 static void wakeup(int sig)
  52 {
  53   if (fcntl(0, F_SETFL, stdin_fdflags | O_NONBLOCK))
  54     barf("fcntl(stdin)", errno);
  55 }
  56
  57 static int intarg(const char *p, const char *what)
  58 {
  59   int e = errno;
  60   long n;
  61   char *q;
  62
  63   errno = 0;
  64   n = std::strtol(p, &q, 0);
  65   if (*q || errno || n < 0 || n > INT_MAX)
  66     barf("bad int %s `%s'", 0, what, p);
  67   errno = e;
  68   return ((int)n);
  69 }
  70
  71 static NTL::GF2X polyarg(const char *p, const char *what)
  72 {
  73   std::string s{p};
  74
  75   // Oh, this is wretched: NTL reads and writes the nibbles backwards.  Hack
  76   // hack.
  77   if (s.size() < 2 || s[0] != '0' || s[1] != 'x') goto bad;
  78   std::reverse(s.begin() + 2, s.end());
  79
  80   { NTL::GF2X a;
  81     std::istringstream ss{s};
  82     ss >> a;
  83     if (!ss) goto bad;
  84     return (a);
  85   }
  86 bad:
  87   barf("bad poly %s `%s'", 0, what, p);
  88 }
  89
  90 static unsigned long long hash(const unsigned char *p, size_t n)
  91 {
  92   size_t i;
  93   unsigned long long h = 0x6a078c523ae42f68ull;
  94
  95   for (i = 0; i < n; i++) { h += p[i]; h *= 0xbeaa913866b8d5a3ull; }
  96   return (h);
  97 }
  98
  99 static std::string showpoly(NTL::GF2X f)
 100 {
 101   std::ostringstream ss;
 102   ss << f;
 103   std::string s{ss.str()};
 104   assert(s.size() >= 2 && s[0] == '0' && s[1] == 'x');
 105   std::reverse(s.begin() + 2, s.end());
 106   return (s);
 107 }
 108
 109 static NTL::ZZ zzarg(const char *p, const char *what)
 110 {
 111   std::string s{p};
 112   std::istringstream ss{s};
 113   NTL::ZZ n;
 114   ss >> n;
 115   if (!ss) barf("bad int %s `%s'", 0, what, p);
 116   return (n);
 117 }
 118
 119 static void seed()
 120 {
 121   unsigned char b[256];
 122   FILE *fp;
 123
 124   if ((fp = std::fopen("/dev/urandom", "rb")) == 0)
 125     barf("open(/dev/urandom)", errno);
 126   if (std::fread(b, 1, sizeof(b), fp) != sizeof(b)) {
 127     barf("read(/dev/urandom)%s",
 128          ferror(fp) ? errno : 0,
 129          ferror(fp) ? "" : ": unexpected short read");
 130   }
 131   std::fclose(fp);
 132   NTL::ZZ t{NTL::ZZFromBytes(b, sizeof(b))};
 133   SetSeed(t);
 134 }
 135
 136 struct step {
 137   NTL::ZZ_p du, dv;
 138   NTL::GF2E m;
 139 };
 140
 141 struct hist {
 142   unsigned long long k;
 143   NTL::ZZ_p u, v;
 144   NTL::GF2E y;
 145 };
 146
 147 #define NSTEP 23
 148 #define NSQR 2
 149 #define NHIST 8
 150
 151 #define CHECK_NITER 65536
 152
 153 int main(int argc, char *argv[])
 154 {
 155   int i;
 156   int dpbits;
 157   unsigned char buf[4096];
 158   unsigned long long niter, dpmask;
 159   fd_set fds_in;
 160   struct timeval tv;
 161   struct sigaction sa;
 162   ssize_t n;
 163
 164   prog = argv[0];
 165
 166   if (argc != 7) {
 167     std::fprintf(stderr, "usage: %s DPBITS gf2x P A B L\n",
 168                  prog);
 169     std::exit(2);
 170   }
 171   dpbits = intarg(argv[1], "dpbits");
 172   dpmask = (1ull << dpbits) - 1;
 173   if (std::strcmp(argv[2], "gf2x") != 0)
 174     barf("unknown group kind `%s'", 0, argv[1]);
 175   NTL::GF2X p = polyarg(argv[3], "p");
 176   if (!IterIrredTest(p)) barf("not irreducible", 0);
 177   NTL::GF2E::init(p);
 178
 179   NTL::GF2E a{to_GF2E(polyarg(argv[4], "a"))};
 180   if (a == 0 || a == 1) barf("a is trivial", 0);
 181   NTL::GF2E b{to_GF2E(polyarg(argv[5], "b"))};
 182   if (b == 0) barf("b isn't a unit", 0);
 183   NTL::ZZ l{zzarg(argv[6], "l")};
 184   if (!ProbPrime(l)) barf("order isn't prime", 0);
 185   NTL::ZZ_p::init(l);
 186
 187   NTL::GF2X::HexOutput = 1;
 188
 189   if (power(a, l) != 1) barf("a has wrong order", 0);
 190   if (power(b, l) != 1) barf("b has wrong order", 0);
 191
 192   // Build the table of steps.  Must do this deterministically!
 193   step tab[NSTEP - NSQR];
 194   SetSeed(NTL::ZZ::zero());
 195   for (i = 0; i < NSTEP - NSQR; i++) {
 196     tab[i].du = NTL::random_ZZ_p();
 197     tab[i].dv = NTL::random_ZZ_p();
 198     tab[i].m = power(a, rep(tab[i].du))*power(b, rep(tab[i].dv));
 199   }
 200
 201   // Now we start the random walk...
 202   seed();
 203   niter = 8ull << (dpbits ? dpbits : (NumBits(l) + 1)/2);
 204 again:
 205   NTL::ZZ_p u{NTL::random_ZZ_p()}, v{NTL::random_ZZ_p()};
 206   NTL::GF2E t = power(a, rep(u))*power(b, rep(v));
 207
 208   hist h[NHIST];
 209   unsigned o = 0;
 210   unsigned long long k = 0;
 211
 212   if (!dpbits)
 213     for (i = 0; i < NHIST; i++) h[i].k = 0;
 214
 215   stdin_fdflags = fcntl(0, F_GETFL);
 216   if (stdin_fdflags < 0) barf("fcntl stdin", errno);
 217   sa.sa_handler = cleanup;
 218   sigemptyset(&sa.sa_mask);
 219   sigaddset(&sa.sa_mask, SIGTSTP);
 220   sigaddset(&sa.sa_mask, SIGCONT);
 221   sa.sa_flags = 0;
 222   if (sigaction(SIGINT, &sa, 0) ||
 223       sigaction(SIGTERM, &sa, 0) ||
 224       sigaction(SIGQUIT, &sa, 0) ||
 225       sigaction(SIGTSTP, &sa, 0))
 226     barf("sigaction", errno);
 227   sa.sa_handler = wakeup;
 228   if (sigaction(SIGCONT, &sa, 0))
 229     barf("sigaction", errno);
 230   if (fcntl(0, F_SETFL, stdin_fdflags | O_NONBLOCK))
 231     barf("fcntl stdin", errno);
 232
 233   for (;;) {
 234     if (k >= niter) goto again;
 235     if (!(k%CHECK_NITER)) {
 236       FD_ZERO(&fds_in); FD_SET(0, &fds_in);
 237       tv.tv_sec = 0; tv.tv_usec = 0;
 238       if (select(1, &fds_in, 0, 0, &tv) < 0)
 239         { if (errno != EINTR) barf("select", errno); }
 240       else if (FD_ISSET(0, &fds_in)) {
 241         for (;;) {
 242           n = read(0, buf, sizeof(buf));
 243           if (!n) { cleanup(0); std::exit(1); }
 244           else if (n > 0) continue;
 245           else if (errno == EAGAIN) break;
 246           else if (errno != EINTR) barf("read stdin", errno);
 247         }
 248       }
 249     }
 250     k++;
 251
 252     size_t nb = NumBytes(rep(t));
 253     BytesFromGF2X(buf, rep(t), nb);
 254     unsigned long long hh = hash(buf, nb);
 255
 256     if (dpbits) {
 257       if (!(hh&dpmask)) {
 258         std::cout << u << " " << v << " " << showpoly(rep(t)) << " "
 259                   << k << std::endl;
 260         goto done;
 261       }
 262     } else {
 263       for (i = 0; i < NHIST; i++) {
 264         if (t == h[i].y) {
 265           if (h[i].u == u || h[i].v == v) goto again;
 266           std::cout << (h[i].u - u)/(v - h[i].v) << " " << k << std::endl;
 267           goto done;
 268         }
 269       }
 270
 271       if (k > 3*h[o].k) {
 272         h[o].u = u; h[o].v = v; h[o].y = t; h[o].k = k;
 273         o = (o + 1)%NHIST;
 274       }
 275     }
 276
 277     i = hh%NSTEP;
 278     if (i >= NSTEP - NSQR) { mul(u, u, 2); mul(v, v, 2); sqr(t, t);; }
 279     else { add(u, u, tab[i].du); add(v, v, tab[i].dv); mul(t, t, tab[i].m); }
 280   }
 281
 282 done:
 283   cleanup(0);
 284   return (0);
 285 }