euler-0288.cpp

// ////////////////////////////////////////////////////////
// # Title
// An enormous factorial
//
// # URL
// https://projecteuler.net/problem=288
// http://euler.stephan-brumme.com/288/
//
// # Problem
// For any prime `p` the number `N(p,q)` is defined by `N(p,q) = sum_{n=0 ... q}{T_n * p^n}`
// with `T_n` generated by the following random number generator:
//
// `S_0 = 290797`
// `S_{n+1} = S_n^2 mod 50515093`
// `T_n = S_n mod p`
//
// Let `Nfac(p,q)` be the factorial of `N(p,q)`.
// Let `NF(p,q)` be the number of factors `p` in `Nfac(p,q)`.
//
// You are given that `NF(3,10000) mod 3^20 = 624955285`.
//
// Find `NF(61, 10^7) mod 61^10`
//
// # Solved by
// Stephan Brumme
// October 2017
//
// # Algorithm
// I started with the easy part:
// The function ''countFactors()'' returns the exponent of a given prime in a factorial's prime factorization.
// It relies on Legendre's formula which I already used for other problems, too ( see https://en.wikipedia.org/wiki/Legendre%27s_formula ).
//
// The pseudo-random number generator wasn't hard as well. I wrote that stuff when I saw the problem for the first time but then I was stuck and continued with other problems.
// Upen revisiting this problem I made the crucial observation that the result has to be printed `mod 61^10` ==> and 61 is the same prime used for `NF(61, 10^7)`.
// (the same hold true for the example `NF(3, 10000) mod 3^20`).
//
// Therefore the result depends on three parameters:
// - a prime `p = 61`
// - a positive number `q = 10^7`
// - an exponent of the modulus `exponent = 10`
//
// The term `p^n` in the formula `N(p,q) = sum_{n=0 ... q}{T_n * p^n}` will be `p=61` and `n=0,1,...` therefore `61^0 = 1`, `61^1 = 61`, `61^2 = 3721`, ...
// It can't exceed `61^10` because of the `mod 61^10` part in the final formula. The variable ''maxPower'' will be 1, 61, 3721, ... until it reaches `61^10` and remains at the level.
//
// There are only 61 possible values for `T_n = S_n mod 61`. That means that the number of different `T_n * p^n` is very small (in fact it's 71).
// Caching those results of ''countFactors(t * p^n)'' let's the runtime of the program drop from 2.6 to 0.35 seconds.
//
// # Note
// `T_n` is at most 60, so `max(T_n * 61^10) = 60 * 61^10 approx 4.28 * 10^19` is the largest value I will encounter in my computations.
// Unfortunately this value is beyond the range of 64 bit integers (it requires 66 bits) and I needed to switch to GCC's 128 bit extension.
// Therefore this code doesn't compile with Visual C++ etc.
//
// I didn't like this problem because I spent too much time "discovering" that I need 128 bit arithmetic.

#include <iostream>
#include <vector>
#include <unordered_map>

// I had overflows with 64 bit arithmetic, so I switched to GCC's 128 bit extension
typedef __int128 Number;

// compute the exponent of prime in the factorization of f!
unsigned long long countFactors(Number f, unsigned int prime)
{
  // Legendre's formula, https://en.wikipedia.org/wiki/Legendre%27s_formula
  unsigned long long result = 0;
  Number power  = prime;
  while (power <= f)
  {
    result += f / power;
    power  *= prime;
  }

  return result;
}

// variables of the problem statement:
// p = prime      = 61
// q = iterations = 10^7
// moduloExponent = 10 (because of 61^10)
unsigned long long solve(unsigned int prime, unsigned int iterations, unsigned int exponent)
{
  // compute prime^moduloExponent => 61^10
  unsigned long long modulo = 1;
  for (unsigned int i = 1; i <= exponent; i++)
    modulo *= prime;

  // will be 61^i => 61^0, 61^1, 61^2, ..., 61^10 and keep it at 61^10
  Number maxPower = 1;

  // seed of pseudo-random number generator
  unsigned long long s = 290797;

  unsigned long long result = 0;
  for (unsigned int i = 0; i <= iterations; i++)
  {
    // T(i, 61) * 61^i
    auto t = s % prime;
    // => but 61^i limited to 61^10
    auto product = t * maxPower;

    // I observed only 71 different products (in 1000000 iterations !), let's cache them
    static std::unordered_map<Number, unsigned long long> cache;
    auto lookup = cache.find(product);
    if (lookup == cache.end())
    {
      // new, unknown parameters, must call countFactors()
      auto current   = countFactors(product, prime);
      // cache the value
      cache[product] = current;
      // add it to the result, too
      result        += current;
    }
    else
    {
      // add looked up value
      result += lookup->second;
    }

    // keep under 61^10
    result %= modulo;

    // 61^i => 61^(i+1), but limit to 61^10
    if (maxPower < modulo)
      maxPower *= prime;

    // next random number
    s *= s;
    s %= 50515093;
  }

  return result;
}

int main()
{
  unsigned int prime      = 61;
  unsigned int iterations = 10000000;
  unsigned int exponent   = 10;
  std::cin >> prime >> iterations >> exponent;
  std::cout << solve(prime, iterations, exponent) << std::endl;
  return 0;
}