Computational "proof" for solution of "whiteout" from Tokyo Westerns/MMA CTF 2016

whiteout.cpp

// We want to prove that 5618427494400 is the maximum sigma in range 1..10^12
//
// The following algorithm will explore the entire search space and prune it
// using an upper bound for sigma, while considering the lower bound given by
// our example of sigma(995886571680) = 5618427494400
//
// We consider only primes up to sqrt(10^12). A single prime p is obviously
// not the solution because the sum of its divisors is
// p + 1 < 10^12 < 5618427494400

#include <bits/stdc++.h>
using namespace std;

using ull = unsigned long long;

ull target = 1000000000000ULL;

vector<int> primes, expo;

//ull best_sigma = 1, best_prod = 1;
ull best_sigma = 5618427494400, best_prod = 995886571680;

const ull inf = numeric_limits<ull>::max();

// Multiply and cap at 2^64 - 1
ull multiply(ull a, ull b) {
  return (a <= inf / b) ? a * b : inf;
}

// Compute floor(log_{base}(x))
int log_floor(ull base, ull x) {
  ull prod = 1;
  int res = 0;
  while (prod <= x) {
    prod *= base;
    ++res;
  }
  return res - 1;
}

// Computes an upper bound on sigma(n) for n <= x and n only has prime
// factors >= p.
//
// For a single prime p, we know that sigma(p^e) = (p^(e+1)-1)/(p-1), so
// sigma(p^e) / p^e <= 1 + 1 / (p - 1)
//
// Also for each q > p we have sigma(q^e) / q^e <= 1 + 1 / (p - 1)
//
// Hence sigma(n) / n <= (1 + 1 / (p - 1))^m where m is the number of distinct
// prime factors of n. We have m <= log_p(x), so in total:
//
// sigma(n) <= x * (1 + 1 / (p - 1))^log_p(x)
ull compute_upper(ull x, int p) {
  for (int i = log_floor(p, x); i >= 1; --i)
    x = (multiply(x + 1, p) - 2) / (p - 1);
  return x;
}

// Backtrack to the optimal solution
//
// i = index of next prime
// prod = product of chosen prime powers amongst primes with indexes 0..{i-1}
// cur = sigma(prod)
void dfs(size_t i=0, ull prod=1, ull cur=1) {
  if (i == primes.size()) return;
  if (prod > target) return;

  int p = primes[i];

  ull upper = multiply(cur, compute_upper(target / prod, primes[i]));
  if (best_sigma > upper)
    return;

  if (cur > best_sigma) {
    best_sigma = cur;
    best_prod = prod;
  }

  ull x = p;
  for (expo[i] = 0; prod <= target; ++expo[i]) {
    dfs(i + 1, prod, cur * (x - 1) / (p - 1));
    prod *= p;
    x *= p;
  }
}

// Compute primes up to sqrt(target)
void sieve() {
  ull upper = 1;
  while (upper * upper <= target)
    ++upper;
  vector<bool> sieve(upper + 1);
  for (ull i = 2; i <= upper; ++i) {
    if (!sieve[i])
      primes.push_back(i);
    for (ull j = i*i; j <= upper; j += i)
      sieve[j] = 1;
  }
}

int main() {
  sieve();
  expo.resize(primes.size());
  dfs();
  cout << best_prod << " " << best_sigma << endl;
}