roll.cpp (ver9)

/*
  This is a port of the roll.py Python 3 library to C++.

  It works as long as the individual numbers of the input are less than 2^32.

  The main function computes the divisors less than 10^8 of 10 tri 10 + 23.

  See roll.py for more documentation.
*/

#include <cstdint>
#include <vector>
#include <map>
#include <iostream>

using namespace std;

typedef uint64_t nat;

struct roll {
  nat flat;
  nat loop;
};

nat operator%(nat a, roll b) {
  return a < b.flat ? a : (a - b.flat) % b.loop + b.flat;
}

struct divisor_cache {
  vector<nat> primes;
  vector<nat> phi;
};

divisor_cache make_cache(nat n) {
  vector<bool> old(n);  // sieve of eratosthenes
  vector<nat> primes;
  vector<nat> phi(n, 1);
  phi[0] = 0;
  for (size_t i = 2; i < n; i++) {
    if (not old[i]) {
      primes.push_back(i);
      for (size_t j = i; j < n; j += i) {
        old[j] = true;
        phi[j] *= i-1;
      }
      for (nat k = i*i; k < n; k *= i) {
        for (size_t j = k; j < n; j += k) {
          phi[j] *= i;
        }
      }
    }
  }
  return {primes, phi};
}

nat pow_mod(nat a, nat b, roll c) {
  nat d = 1;
  while (b) {
    if (b % 2) d = d * a % c;
    a = a * a % c;
    b /= 2;
  }
  return d;
}

nat gcd(nat a, nat b) {
  if (a == 0 and b == 0) {
    return 0;
  }
  while (b) {
    nat c = a % b;
    a = b;
    b = c;
  }
  return a;
}

roll roll_of_powers(nat a, roll b, const divisor_cache& cache) {
  if (a == 0) return {1, 1};
  if (a == 1) return {0, 1};
  nat flat = 0;
  nat a_pow_flat = 1;
  while (a_pow_flat < b.flat) {
    flat++;
    a_pow_flat *= a;
  }
  nat d = b.loop / gcd(b.loop, a_pow_flat);
  nat c = gcd(d, a);
  while (c != 1) {
    d /= c;
    c = gcd(d, c);
    flat++;
  }
  return {flat, cache.phi[d]};
}

// (tower of exponents a) % (roll b)
nat mod(const vector<nat>& a, roll b, const divisor_cache& cache) {
  if (a.size() == 0) return 1 % b;
  vector<roll> rolls(a.size());
  rolls[0] = b;
  for (size_t i = 0; i < a.size()-1; i++) {
    rolls[i+1] = roll_of_powers(a[i], rolls[i], cache);
  }
  nat c = 1 % b;
  for (size_t i = a.size()-1; i != (size_t) -1; i--) {
    c = pow_mod(a[i], c, rolls[i]);
  }
  return c;
}

void print_divisors(const vector<nat>& tower, nat term, const divisor_cache& cache) {
  for (nat p : cache.primes) {
    if ((mod(tower, {0, p}, cache) + term) % p == 0) {
      cout << p << " " << flush;
    }
  }
}

int main(int argc, char ** argv) {
  nat n = 100000000;
  divisor_cache cache = make_cache(n);
  vector<nat> tower(10, 10);
  print_divisors(tower, 23, cache);
  cout << endl;
}