#include <complex>
#include <fftw3.h>
#include <github.com/apronchenkov/perfmon/public/perfmon.h>
#include <iostream>
#include <vector>

using Real = double;
using Complex = std::complex<Real>;
using FFTWComplex = fftw_complex;
using FFTWPlan = fftw_plan;
template <class... Args>
auto FFTWSetTimelimit(Args&&... args)
    -> decltype(fftw_set_timelimit(std::forward<Args>(args)...)) {
  return fftw_set_timelimit(std::forward<Args>(args)...);
}
template <class... Args>
auto FFTWPlanDFT1D(Args&&... args)
    -> decltype(fftw_plan_dft_1d(std::forward<Args>(args)...)) {
  return fftw_plan_dft_1d(std::forward<Args>(args)...);
}
template <class... Args>
auto FFTWDestroyPlan(Args&&... args)
    -> decltype(fftw_destroy_plan(std::forward<Args>(args)...)) {
  return fftw_destroy_plan(std::forward<Args>(args)...);
}
template <class... Args>
auto FFTWExecutePlan(Args&&... args)
    -> decltype(fftw_execute(std::forward<Args>(args)...)) {
  return fftw_execute(std::forward<Args>(args)...);
}
//*/
/*
using Real = float;
using Complex = std::complex<Real>;
using FFTWComplex = fftwf_complex;
using FFTWPlan = fftwf_plan;
template <class... Args>
auto FFTWSetTimelimit(Args&&... args)
    -> decltype(fftw_set_timelimit(std::forward<Args>(args)...)) {
  return fftwf_set_timelimit(std::forward<Args>(args)...);
}
template <class... Args>
auto FFTWPlanDFT1D(Args&&... args)
    -> decltype(fftwf_plan_dft_1d(std::forward<Args>(args)...)) {
  return fftwf_plan_dft_1d(std::forward<Args>(args)...);
}
template <class... Args>
auto FFTWDestroyPlan(Args&&... args)
    -> decltype(fftwf_destroy_plan(std::forward<Args>(args)...)) {
  return fftwf_destroy_plan(std::forward<Args>(args)...);
}
template <class... Args>
auto FFTWExecutePlan(Args&&... args)
    -> decltype(fftwf_execute(std::forward<Args>(args)...)) {
  return fftwf_execute(std::forward<Args>(args)...);
}
//*/

class FFTEng {
 public:
  FFTEng() = default;

  explicit FFTEng(int n) : FFTEng(n, FFTW_ESTIMATE) {}

  FFTEng(int n, unsigned int fftwFlag) {
    PERFMON_SCOPE("FFTEng.ctor");
    x_.resize(n);
    y_.resize(n);
    forward_ = FFTWPlanDFT1D(n, reinterpret_cast<FFTWComplex*>(x_.data()),
                             reinterpret_cast<FFTWComplex*>(y_.data()),
                             FFTW_FORWARD, fftwFlag);
    backward_ = FFTWPlanDFT1D(n, reinterpret_cast<FFTWComplex*>(y_.data()),
                              reinterpret_cast<FFTWComplex*>(x_.data()),
                              FFTW_BACKWARD, fftwFlag);
    if (!forward_ || !backward_) {
      if (forward_) {
        FFTWDestroyPlan(forward_);
        forward_ = nullptr;
      }
      if (backward_) {
        FFTWDestroyPlan(backward_);
        backward_ = nullptr;
      }
      throw std::runtime_error("fftw_plan_dft_1d: Failed.");
    }
  }

  FFTEng(FFTEng&& rhs) : FFTEng() {
    std::swap(x_, rhs.x_);
    std::swap(y_, rhs.y_);
    std::swap(forward_, rhs.forward_);
    std::swap(backward_, rhs.backward_);
  }

  FFTEng(const FFTEng& rhs) = delete;

  FFTEng& operator=(FFTEng&& rhs) {
    std::swap(x_, rhs.x_);
    std::swap(y_, rhs.y_);
    std::swap(forward_, rhs.forward_);
    std::swap(backward_, rhs.backward_);
    return *this;
  }

  FFTEng& operator=(const FFTEng& rhs) = delete;

  ~FFTEng() {
    if (forward_) {
      FFTWDestroyPlan(forward_);
      forward_ = nullptr;
    }
    if (backward_) {
      FFTWDestroyPlan(backward_);
      backward_ = nullptr;
    }
  }

  void Forward() {
    PERFMON_SCOPE("FFTEng.Forward");
    FFTWExecutePlan(forward_);
  }

  void Backward() {
    PERFMON_SCOPE("FFTEng.Backward");
    FFTWExecutePlan(backward_);
    const Real scale = 1.0 / x_.size();
    for (auto& x : x_) {
      x *= scale;
    }
  }

  int size() const { return x_.size(); }

  Complex& X(size_t i) { return x_[i]; }
  const Complex& X(size_t i) const { return x_[i]; }
  auto XBegin() { return x_.begin(); }
  auto XBegin() const { return x_.begin(); }
  auto XEnd() { return x_.end(); }
  auto XEend() const { return x_.end(); }

  Complex& Y(size_t i) { return y_[i]; }
  const Complex& Y(size_t i) const { return y_[i]; }
  auto YBegin() { return y_.begin(); }
  auto YBegin() const { return y_.begin(); }
  auto YEnd() { return y_.end(); }
  auto YEend() const { return y_.end(); }

 private:
  std::vector<Complex> x_;
  std::vector<Complex> y_;
  FFTWPlan forward_{nullptr};
  FFTWPlan backward_{nullptr};
};

std::vector<int> GenPrimes(int nmax) {
  std::vector<int> result;
  std::vector<bool> sieve(nmax, true);
  sieve[0] = false;
  sieve[1] = false;
  for (int i = 2; i <= nmax; ++i) {
    if (!sieve[i]) {
      continue;
    }
    result.push_back(i);
    for (int j = i + i; j <= nmax; j += i) {
      sieve[j] = false;
    }
  }
  return result;
}

template <class Container>
bool Any(const Container& input) {
  for (const auto& x : input) {
    if (x) {
      return true;
    }
  }
  return false;
}

void PrintBits(const std::vector<bool>& input) {
  const size_t kMaxSize = 100;
  for (size_t i = 0; i < kMaxSize && i < input.size(); ++i) {
    std::cout << input[i];
  }
  if (input.size() > kMaxSize) {
    std::cout << "...";
  }
}

int main() {
  const int kNMax = 200000;
  FFTWSetTimelimit(.1);
  FFTEng eng(407680);

  std::fill(eng.XBegin(), eng.XEnd(), 0);
  const auto primes = GenPrimes(kNMax);
  for (size_t i = 0; i < primes.size() && primes[i] <= kNMax; ++i) {
    eng.X(primes[i]) = 1;
    for (size_t j = i; j < primes.size() && primes[j] <= kNMax / primes[i];
         ++j) {
      eng.X(primes[i] * primes[j]) = 1;
    }
  }
  eng.Forward();
  const std::vector<Complex> aY(eng.YBegin(), eng.YEnd());

  std::vector<bool> candidates(kNMax, true);
  for (int n = 0; Any(candidates); ++n) {
    std::vector<bool> isUndefined = candidates;
    std::vector<bool> g(isUndefined.size());
    while (Any(isUndefined)) {
      std::fill(eng.XBegin(), eng.XEnd(), 0);
      for (size_t i = 0; i < isUndefined.size(); ++i) {
        eng.X(i) = Complex(isUndefined[i], g[i]);
      }
      eng.Forward();
      for (size_t i = 0; i < aY.size(); ++i) {
        eng.Y(i) *= aY[i];
      }
      eng.Backward();
      for (size_t i = 0; i < isUndefined.size(); ++i) {
        isUndefined[i] = candidates[i] && (eng.X(i).real() > 0.5) &&
                         !(eng.X(i).imag() > 0.5);
        g[i] = candidates[i] && !(eng.X(i).real() > 0.5) &&
               !(eng.X(i).imag() > 0.5);
      }
    }
    std::cout << n << ":\t";
    PrintBits(g);
    std::cout << '\n';

    std::fill(eng.XBegin(), eng.XEnd(), 0);
    for (size_t i = 0; i < g.size(); ++i) {
      eng.X(i) = Complex(0, g[i]);
    }
    eng.Forward();
    for (size_t i = 0; i < aY.size(); ++i) {
      eng.Y(i) *= aY[i];
    }
    eng.Backward();
    for (int i = 0; i < kNMax; ++i) {
      candidates[i] = candidates[i] && !g[i] && (eng.X(i).imag() > 0.5);
    }
  }

  for (const auto& counter : PERFMON_COUNTERS()) {
    std::cerr << counter.Name() << ": " << counter.Calls() << ' '
              << counter.AverageSeconds() << '\n';
  }

  return 0;
}