supercomputer-solver.cc

#include <cassert>

#include <NTL/mat_GF2.h>
#include <NTL/GF2.h>

using NTL::GF2;
using NTL::mat_GF2;
using NTL::conv;

int g(GF2 x1, GF2 z1, GF2 x2, GF2 z2) {
  size_t x = 0;
  if (IsOne(x1))  x += 2;
  if (IsOne(z1))  x += 1;
  size_t y = 0;
  if (IsOne(x2))  y += 2;
  if (IsOne(z2))  y += 1;
  int table[4][4] = {
    {0, 0, 0, 0},
    {0, 0, 1, -1},
    {0, -1, 0, 1},
    {0, 1, -1, 0}
  };
  return table[x][y];
}

struct MeasurementResult {
  GF2 result;
  bool was_random;
};

class CircuitRunner {
  public:
    CircuitRunner(size_t bits) :
      state_(NTL::ident_mat_GF2(2*bits + 1)),
      bits_(bits) { }
    void Hadamard(size_t a) {
      for (size_t i = 0; i < 2*bits_; i++) {
        auto& row = state_[i];
        row[2*bits_] += row[a] * row[a+bits_];
        swap(row[a], row[a+bits_]);
      }
    }
    void Phase(size_t a) {
      for (size_t i = 0; i < 2*bits_; i++) {
        auto& row = state_[i];
        row[2*bits_] += row[a] * row[a+bits_];
        row[a+bits_] += row[a];
      }
    }
    void CNot(size_t a, size_t b) {
      // a is control, b is target
      for (size_t i = 0; i < 2*bits_; i++) {
        auto& row = state_[i];
        row[2*bits_] += row[a] * row[b+bits_] *
          (row[a+bits_] + row[b] + conv<GF2>(1));
        row[b] += row[a];
        row[a+bits_] += row[b+bits_];
      }
    }
    MeasurementResult Measure(size_t a) {
      for (size_t p = bits_; p < 2*bits_; p++) {
        if (IsZero(state_[p][a]))  continue;
        /* Random case */
        for (size_t i = 0; i < 2*bits_; i++) {
          if (i == p)  continue;
          if (IsZero(state_[i][a]))  continue;
          RowSum(i, p);
        }
        state_[p-bits_] = state_[p];
        NTL::clear(state_[p]);
        GF2 rand = conv<GF2>(1);
        state_[p][2*bits_] = rand;
        state_[p][a+bits_] = conv<GF2>(1);
        return {rand, true};
      }
      /* Deterministic case */
      NTL::clear(state_[2*bits_]);
      for (size_t i = 0; i < bits_; i++) {
        if (IsZero(state_[i][a]))  continue;
        RowSum(2*bits_, i+bits_);
      }
      return {state_[2*bits_][2*bits_], false};
    }
    void AssertIs(size_t a, GF2 result) {
      MeasurementResult m = Measure(a);
      assert(!m.was_random);
      assert(m.result == result);
    }
    void AssertZero(size_t a) {
      MeasurementResult m = Measure(a);
      assert(!m.was_random);
      assert(IsZero(m.result));
    }
    void AssertOne(size_t a) {
      MeasurementResult m = Measure(a);
      assert(!m.was_random);
      assert(!IsZero(m.result));
    }
    GF2 MeasureAssertRandom(size_t a) {
      MeasurementResult m = Measure(a);
      assert(m.was_random);
      return m.result;
    }
    GF2 MeasureCheckDeterministic(size_t a) {
      MeasurementResult m = Measure(a);
      if (m.was_random) {
        printf("%lu was random\n", a);
      }
      return m.result;
    }


  private:
    void RowSum(size_t h, size_t i) {
      // onto h, from i
      int accum = 0;
      if (IsOne(state_[h][2*bits_])) {
        accum += 2;
      }
      if (IsOne(state_[i][2*bits_])) {
        accum += 2;
      }
      for (size_t j = 0; j < bits_; j++) {
        accum += g(
            state_[i][j], state_[i][j+bits_],
            state_[h][j], state_[h][j+bits_]);
      }
      GF2 rh = conv<GF2>((accum%4 == 0) ? 0 : 1);

      state_[h] += state_[i];
      state_[h][2*bits_] = rh;
    }
    mat_GF2 state_;
    size_t bits_;
};

void zero_test(void) {
  CircuitRunner c(3);
  c.AssertZero(0);
  c.AssertZero(0);
  c.AssertZero(1);
  c.AssertZero(0);
  c.AssertZero(2);
  c.AssertZero(1);
  c.AssertZero(0);
}


void phase_test(void) {
  CircuitRunner c(3);
  c.Hadamard(0);
  c.CNot(0, 1);
  c.Phase(0);
  c.CNot(0, 2);
  c.CNot(0, 1);
  c.Phase(0);
  c.CNot(0, 2);
  c.Hadamard(0);
  c.AssertOne(0);
  c.AssertZero(1);
  c.AssertZero(2);
}

void bellpair_test(void) {
  CircuitRunner c(2);
  c.Hadamard(0);
  c.CNot(0, 1);
  c.AssertIs(1, c.MeasureAssertRandom(0));
}

uint64_t doit(void) {
  CircuitRunner c(48);
#include "circuit.inl"
  uint64_t result = 0;
  for (int i = 0; i < 48; i++) {
    result *= 2;
    if (!IsZero(c.Measure(i).result)) {
      result++;
    }
  }
  return result;
}


int main(void) {
  zero_test();
  hadamard_test();
  phase_test();
  bellpair_test();
  printf("FLAG{%lu}\n", doit());
  return 0;
}