euler-0375.cpp

// ////////////////////////////////////////////////////////
// # Title
// Minimum of subsequences
//
// # URL
// https://projecteuler.net/problem=375
// http://euler.stephan-brumme.com/375/
//
// # Problem
// Let `S_n` be an integer sequence produced with the following pseudo-random number generator:
// `S_0     = 290797`
// `S_{n+1} = S^2_n mod 50515093`
//
// Let `A(i, j)` be the minimum of the numbers `S_i`, `S_{i+1}`, ... , `S_j` for `i <= j`.
// Let `M(N) = sum{A(i, j)}` for `1 <= i <= j <= N`.
// We can verify that `M(10) = 432256955` and `M(10 000) = 3264567774119`.
//
// Find `M(2 000 000 000)`.
//
// # Solved by
// Stephan Brumme
// December 2017
//
// # Algorithm
// I wrote a ''bruteForce()'' algorithm in a few minutes which - to my surprise - even solves the `M(10000)` case in less than a second.
// After precomputing the first 10000 pseudo-random values it performs the following tasks for each position:
// - set ''minimum'' to the current value ''data[i]''
// - go backwards until reaching the first position, compare ''minimum'' to ''data[j]'' and update it according
// - add ''minimum'' to the result
//
// This `O(n^2)` algorithm repeatedly processes the same values.
// It's more efficient to keep track of the positions where ''data[j]'' was smaller than ''minimum''.
// For example, if the next ''minimum'' is 5 positions away, then I can easily add 5 times the current ''minimum'' to ''result''.
//
// The ''search()'' function implements this idea:
// - it has a stack ''best'' which contains those numbers smaller than the current number ("previous ''minimum's") and their position in the stream of pseudo-random numbers
// - this stack is initialized with ''0'', a number smaller than anything generated by ''pseudoRandom'' so that the stack is never empty (prevents a few corner-cases)
// - whenever a new pseudo-random number called ''current'' is generated then I reduce the stack until the top-most element is smaller than ''current''
// - ''current'' and its position ''i'' are pushed to the stack
// - for all elements of the stack I add their value and their distance to the next element of the stack to ''result''
//
// # Note
// The stack ''best'' remains quite small. It contains less than 40 values at any time.
//
// # Alternative
// My program needs just under a minute to solve the problem.
// I didn't realize that ''pseudoRandom()'' has a short cycle:
// - there can be at most 50515093 different outputs (50515093 is the modulo)
// - each number depends only on the previous number
// - whenever the internal state becomes a number I have already seen before then a new cycle starts
// - if you enable ''FIND_CYCLE'' then my code will almost instantly display 6308947.
//
// Pretty much everyone on the Project Euler forum spotted the cycle - I didn't ...
// ... but it can dramatically reduce execution time:
// - process the first cycle the same way I did
// - remember the value of ''result''
// - process the second cycle and determine how much ''result'' changed (let's call it ''delta'')
// - skip 6308947 iterations of the loop and just add ''delta'' to the ''result''
// - the last cycle will not be complete, process it the same way I did
//
// You will process less than `2 * 10^6` pseudo-random values instead of `2 * 10^9` and will find the result in less than a second (remember, my code takes a minute).
// A major drawback is that much more code is required for this faster algorithm.

#include <iostream>
#include <vector>
#include <algorithm>

// produce S1, S2, etc.
unsigned int pseudoRandom()
{
  static unsigned long long state = 290797;
  state *= state;
  state %= 50515093;
  return (unsigned int)state;
}

// basic O(n^2) algorithm, surprisingly fast for size = 10000
unsigned long long bruteForce(unsigned int size)
{
  // precompute all values
  std::vector<unsigned int> data = { 0 }; // dummy element to have S1,S2,... at positions 1,2,...
  for (unsigned int i = 1; i <= size; i++)
    data.push_back(pseudoRandom());

  // just a dumb scan over all intervals
  unsigned long long result = 0;
  for (unsigned int i = 1; i <= size; i++)
  {
    auto minimum = data[i];
    //for (unsigned int j = i; j <= size; j++)
    for (unsigned int j = i; j >= 1; j--) // same concept as above but going backwards
    {
      minimum = std::min(minimum, data[j]);
      result += minimum;
    }
  }

  return result;
}

// "number" is the smallest number from its "position" to the current position
struct Minimum
{
  unsigned int number;
  unsigned int position;
};

// almost O(n), about a minute for size = 2*10^9
unsigned long long search(unsigned int size)
{
  // basically a stack of all previous minimums
  std::vector<Minimum> best;
  // add an unused dummy element so that the stack never becomes empty
  best.push_back({ 0, 0 });

  unsigned long long result = 0;
  for (unsigned int i = 1; i <= size; i++)
  {
    // produce next pseudo-random number
    auto current = pseudoRandom();

    // remove all elements larger than the current one
    auto validUntil = i;
    while (best.back().number >= current)
    {
      validUntil = best.back().position;
      best.pop_back();
    }

    // the current value is at least the smallest value for a sequence of length 1 starting at i
    best.push_back({ current, validUntil });

    // add all sequences
    auto last = i + 1;
    for (unsigned int backwards = best.size() - 1; backwards > 0; backwards--)
    {
      // its the same minimum number for all those positions
      unsigned long long distance = last - best[backwards].position;
      result += distance * best[backwards].number;

      // update position
      last = best[backwards].position;
    }
  }

  return result;
}

int main()
{
//#define FIND_CYCLE
#ifdef  FIND_CYCLE
  auto first = pseudoRandom();
  unsigned int cycleLength = 0;
  while (pseudoRandom() != first)
    cycleLength++;
  std::cout << "cycle length: " << cycleLength << std::endl;
#else
  unsigned int limit = 2000000000;
  std::cin >> limit;
  //std::cout << bruteForce(limit) << std::endl;
  std::cout << search(limit) << std::endl;
#endif

  return 0;
}