// ////////////////////////////////////////////////////////
// # 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
#include
#include
// 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 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 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;
}