// ////////////////////////////////////////////////////////
// # Title
// Unfair wager
//
// # URL
// https://projecteuler.net/problem=436
// http://euler.stephan-brumme.com/436/
//
// # Problem
// Julie proposes the following wager to her sister Louise.
// She suggests they play a game of chance to determine who will wash the dishes.
// For this game, they shall use a generator of independent random numbers uniformly distributed between 0 and 1.
// The game starts with `S = 0`.
// The first player, Louise, adds to S different random numbers from the generator until `S > 1` and records her last random number `x`.
// The second player, Julie, continues adding to S different random numbers from the generator until `S > 2` and records her last random number `y`.
// The player with the highest number wins and the loser washes the dishes, i.e. if `y > x` the second player wins.
//
// For example, if the first player draws `0.62` and `0.44`, the first player turn ends since `0.62+0.44 > 1` and `x = 0.44`.
// If the second players draws `0.1`, `0.27` and `0.91`, the second player turn ends since `0.62+0.44+0.1+0.27+0.91 > 2` and `y = 0.91`. Since `y > x`, the second player wins.
//
// Louise thinks about it for a second, and objects: "That's not fair".
// What is the probability that the second player wins?
// Give your answer rounded to 10 places behind the decimal point in the form 0.abcdefghij
//
// # Monte-Carlo simulation written by
// Stephan Brumme
// November 2017
//
// # Idea
// I wrote a simple Monte-Carlo simulation which finds the first 4 digits in a few seconds: 0.5276
//
// My current thinking is that I could implement a better numerical approximation of the probability of drawing `x` as the last number (that is, some `p(x)`).
// Julie's number would be `p(0)` and Louise's `p(p(0))`. The result is the probability that `p(0) < p(p(0))`.
//
// However, I'm not sure whether that would be enough to compute the first 10 digits reliably.
#include
// a simple pseudo-random number generator, returning a value between 0 and 1 (exclusive)
// (produces the same result no matter what compiler you have - unlike rand() from math.h)
double myrand()
{
// based on code from problem 227, modified to return a value between 0 and 1 (exclusive)
static unsigned long long seed = 0;
seed = 6364136223846793005ULL * seed + 1;
const unsigned int ValidBits = 31; // maybe not all bits are random, but it's not a high-quality algorithm anyway
const unsigned int Scaling = 1 << ValidBits;
auto result = ((seed >> 30) % Scaling) / double(Scaling);
return result;
}
double monteCarlo(unsigned int iterations)
{
// count how often player 2 wins
unsigned int twoWins = 0;
for (unsigned int i = 0; i < iterations; i++)
{
// sum of all random numbers
double total = 0;
// last random number of player 1
double lastOne = 0;
while (total < 1)
{
lastOne = myrand();
total += lastOne;
}
// last random number of player 1
double lastTwo = 0;
while (total < 2)
{
lastTwo = myrand();
total += lastTwo;
}
if (lastTwo > lastOne)
twoWins++;
}
return twoWins / double(iterations);
}
int main()
{
// the first 4 digits should be correct - but I have to find the first 10 digits, which is impossible with Monte-Carlo
while (true)
std::cout << monteCarlo(100000000) << std::endl;
return 0;
}