// ////////////////////////////////////////////////////////
// # Title
// Licence plates
//
// # URL
// https://projecteuler.net/problem=371
// http://euler.stephan-brumme.com/371/
//
// # Problem
// Oregon licence plates consist of three letters followed by a three digit number (each digit can be from [0..9]).
// While driving to work Seth plays the following game:
// Whenever the numbers of two licence plates seen on his trip add to 1000 that's a win.
//
// E.g. MIC-012 and HAN-988 is a win and RYU-500 and SET-500 too. (as long as he sees them in the same trip).
//
// Find the expected number of plates he needs to see for a win.
// Give your answer rounded to 8 decimal places behind the decimal point.
//
// Note: We assume that each licence plate seen is equally likely to have any three digit number on it.
//
// # Solved by
// Stephan Brumme
// September 2017
//
// # Algorithm
// The problem description tries hard to confuse the reader ... the plates' letters (e.g. RYU or SET) don't matter at all.
// The only relevant part of a plate is its number. And I've never seen a leading zero on a plate - though it's possible according to the text.
// In the end it doesn't matter as well.
//
// A simple Monte-Carlo simulation can find the first digits (see ''monteCarlo()''):
// - generate a random number ''plate''
// - if ''1000 - plate'' was seen before then the game is over
// - add ''plate'' to the set of seen plates
//
// There are three groups of numbers:
// (a) 500 can only be paired with a second 500
// (b) 0 doesn't pair with any number
// (c) all other numbers
//
// The whole problem can be mapped to states (not the meaning "countries" but "character"):
// - there are 499 pairs in group c ==> 499 states
// - for each of those state there are two options: 500 wasn't seen or it was seen once
// ==> multiplying both results in 998 states
// I decided to create two arrays ''no500'' and ''with500'', each with 499 elements (for a total of 998 values).
//
// And their state transitions:
// - the "zero" plate (group b) doesn't change the current state
//
// The third group (1 - 499 and 501 - 999) has an interesting property:
// for each pair `(x, y)` where `x + y = 1000` is doesn't matter whether `x < y` or `x > y`.
// Therefore I assume that the first number of any pair in group c is always smaller than 500.
// The exact numbers on the plates don't affect the result as long as they belong to group c: observing (12,430,299) is identical to (1,2,3).
// In the end I am only concerned whether a plate's number is new (==> choose the next free number 1,2,3,...), is a duplicate or pairs with an existing number.
// The probabilities are (assuming ''have'' denoted the number of unpaired plates in group c):
// (1) `p_{new} = (998 - 2 * have) / 1000`
// (2) `p_{finished} = have / 1000`
// (3) `p_{duplicate} = have / 1000`
// (4) `p_{zero} = 1 / 1000`
// (5) `p_{500} = 1 / 1000`
//
// Stupid me needs to verify that I got everything right: `1 = p_{new} + p_{finished} + p_{duplicate} + p_{zero} + p_{500}` ==> it's okay !
//
// If I progress from one state to another then the expected value is
// (6) `E = 1 + p_{E_1} * E_1 + p_{E_2} * E_2 + p_{E_3} * E_3 + ...`
//
// I have to be careful with state transitions that lead back to the same state they started from:
// (7) `E = 1 + p_E * E`
// (8) `E - p_E * E = 1`
// (9) `E * (1 - p_E) = 1`
// (10) `E = 1 / (1 - p_E)`
// ==> e.g. assume I throw a fair coin `p_E = 0.5` => `E = 1 / (1 - 0.5) = 1 / 0.5 = 2` so I expect to need two attempts to see the other side
// ==> e.g. assume I throw a biased coin `p_E = 0.6` that has a higher probability to show the last side => `E = 1 / (1 - 0.6) = 1 / 0.4 = 2.5` attempts to see the other side
//
// And the generalized equation is `E = (1 + p_{E_1} * E_1 + p_{E_2} * E_2 + ...) / (1 - p_E)`
//
// I found it much easier to reverse the problem:
// if I have seen all plates 1 to 499 (and there was no pair yet), then
// 1. there is a chance `1 / 1000` to see "zero"
// 2. there is a chance `499 / 1000` to see an old plate
// 3. there is a chance `499 / 1000` to finally see a matching plate
// 4. there is a chance `1 / 1000` to see "500"
//
// The first two cases keep the player in its current state.
// The game is won if case 3 happens-
// If the current state includes that "500" was seen before, then the game is over, too.
// If the current state includes that "500" was __NOT__ seen before, then the state changes to "have500 with 499 plates".
//
// Now the expected number for ''have500[499]'' is:
// (11) `E_{499,have500} = 1 / (1 - (p_{zero} + p_{duplicate})) = 2`
//
// And ''no500[499]'':
// (12) `E_{499,no500} = (1 + 1/1000 * E_{499,have500}) / (1 - (p_{zero} + p_{duplicate})) = 2.4`
//
// Now `E_{498,have500}` can be computed by using the probabilities of staying in the current state and both state transitions (finish or advance to 499 plates).
// With similar argument `E_{498,no500}` depends only on `E_{498,have500}` and `E_{499,no500}`.
// A ''for''-loop computes all values.
//
// Initially, the player hasn't seen a plate at all. That means `E_{0,no500}` is the result.
//
// # Note
// ''search()'' computes bad results if the number of plates is odd because then there is no equivalent to the unique features of the "500" plates.
//
// It took me a few minutes to realize that "zero" is special: first, I didn't think it's a valid plate number at all and second, I'm dumb. Sad but true ...
//
// # Alternative
// A Dynamic Programming approach sounds viable but I highly doubt it can match the efficiency of the plain ''for''-loop.
#include
#include
#include
// a simple pseudo-random number generator
// (produces the same result no matter what compiler you have - unlike rand() from math.h)
unsigned int myrand()
{
// copied from problem 227
static unsigned long long seed = 0;
seed = 6364136223846793005ULL * seed + 1;
return (unsigned int)(seed >> 30);
}
// simulate random journeys ...
double monteCarlo(unsigned int numPlates, unsigned int iterations)
{
unsigned int numPlatesSeen = 0;
for (unsigned int i = 0; i < iterations; i++)
{
// true if plate was seen
std::vector seen(numPlates, false);
while (true)
{
// generate a new plate
auto plate = myrand() % numPlates;
numPlatesSeen++;
// look for matching plate (so that sum is 1000)
if (plate != 0) // 0 isn't matched by any other plate
{
auto other = numPlates - plate;
if (seen[other])
break;
}
// mark current plate as "seen"
seen[plate] = true;
}
}
return numPlatesSeen / double(iterations);
}
// and now smarter, faster, more precise ...
double search(unsigned int numPlates)
{
// how I compute the expected number:
// E = 1 + fail * E
// E - fail * E = 1
// E * (1 - fail) = 1
// E = 1 / (1 - fail)
// if there is another state transition E2 then with similar argument
// E = 1 + fail * E + good * E2
// E = (1 + good * E2) / (1 - fail)
// Seth may observe at most 499 "ordinary" plates
// => zero and 500 are "special"
const auto maxHave = numPlates / 2 - 1;
// several times I have to divide by numPlates and have to make sure the result is a double,
// so let's convert it to double once and keep using plates instead of numPlates from here on
const double plates = numPlates; // actually I thought about invPlates = 1 / double(numPlates)
// track expected number of plates you still have to see if n unpaired plates were already observed
// there is a significant difference whether the "500" plate was seen or not !
std::vector have500(maxHave + 1, 0);
std::vector no500 (maxHave + 1, 0);
// probabilities:
// observe an "old" plate, nothing changes (but the expected number)
auto probDuplicate = maxHave / plates;
// observe a "000" plate, nothing changes (but the expected number)
auto probZero = 1 / plates;
// observe a "500" plate, stop if already saw a "500", otherwise store it
auto prob500 = 1 / plates;
// solve the base case: all 499 plates already seen
// compute probability that the next plate still doesn't finish the game
auto probUnchanged = probDuplicate + probZero;
// now find expected number in case Seth 499 plates plus the "500" plate
have500[maxHave] = 1 / (1 - probUnchanged);
// and find expected number in case Seth 499 plates BUT NOT the "500" plate
no500 [maxHave] = (1 + prob500 * have500[maxHave]) / (1 - probUnchanged);
// continue the computation for 498, 497, ..., 1, 0 observed plates
for (int have = maxHave - 1; have >= 0; have--)
{
// determine number of unpaired plates (exclude zero and 500)
auto numNew = plates - 2*have - 2;
// probability that the plate wasn't seen before but pairs with no other plate yet
auto probNew = numNew / plates;
// it becomes less probable to see an "old" plate ...
probDuplicate = have / plates; // same as probDuplicate -= 1 / plates
probUnchanged = probDuplicate + probZero; // same as probUnchanged += 1 / plates
// same as above: E = (1 + good * E2) / (1 - fail)
have500[have] = (1 + probNew * have500[have + 1]) / (1 - probUnchanged);
no500 [have] = (1 + prob500 * have500[have] + probNew * no500 [have + 1]) / (1 - probUnchanged);
}
// no plate seen yet, not even "500"
return no500[0];
}
int main()
{
unsigned int numPlates = 1000;
std::cin >> numPlates;
std::cout << std::fixed << std::setprecision(8);
std::cout << search(numPlates) << std::endl;
// run simulation
//while (true)
//std::cout << monteCarlo(numPlates, 1000000) << std::endl;
return 0;
}