#include <cmath>
#include <iomanip>
#include <iostream>
#include <random>
#include <utility>
#include <vector>

using namespace std;

double normalCdf(double x)
{
    return 0.5 * erfc(-x / sqrt(2.0));
}

double simulateThresholdStrategy(int a, int b, long long trials, double sigma, mt19937_64& rng)
{
    uniform_int_distribution<int> revealDist(0, 1);
    normal_distribution<double> thresholdDist(0.0, sigma);

    long long wins = 0;
    for (long long i = 0; i < trials; ++i)
    {
        int revealed = (revealDist(rng) == 0) ? a : b;
        double t = thresholdDist(rng);

        bool guessLarger = (revealed > t);
        bool correct = (revealed == b && guessLarger) || (revealed == a && !guessLarger);
        if (correct)
        {
            ++wins;
        }
    }

    return static_cast<double>(wins) / static_cast<double>(trials);
}

double simulateRandomGuess(long long trials, mt19937_64& rng)
{
    uniform_int_distribution<int> coin(0, 1);
    long long wins = 0;
    for (long long i = 0; i < trials; ++i)
    {
        if (coin(rng) == 1)
        {
            ++wins;
        }
    }
    return static_cast<double>(wins) / static_cast<double>(trials);
}

double simulateAlwaysLarger(long long trials, mt19937_64& rng)
{
    uniform_int_distribution<int> revealDist(0, 1);
    long long wins = 0;
    for (long long i = 0; i < trials; ++i)
    {
        if (revealDist(rng) == 1)
        {
            ++wins;
        }
    }
    return static_cast<double>(wins) / static_cast<double>(trials);
}

int main()
{
    const long long trials = 2000000;
    const double sigma = 10.0;

    vector<pair<int, int>> tests = {
        {-3, 3},
        {0, 1},
        {-8, -2},
        {2, 9},
        {-1, 6},
        {1000, 1001}
    };

    random_device rd;
    mt19937_64 rng(rd());

    cout << fixed << setprecision(6);
    cout << "Trials per test = " << trials << "\n";
    cout << "Threshold T ~ N(0, " << sigma * sigma << ")\n\n";

    for (const auto& [a, b] : tests)
    {
        double empirical = simulateThresholdStrategy(a, b, trials, sigma, rng);

        double fa = normalCdf(static_cast<double>(a) / sigma);
        double fb = normalCdf(static_cast<double>(b) / sigma);
        double theoretical = 0.5 + 0.5 * (fb - fa);

        double randomGuess = simulateRandomGuess(trials, rng);
        double alwaysLarger = simulateAlwaysLarger(trials, rng);

        cout << "a=" << a << ", b=" << b << "\n";
        cout << "  threshold strategy (empirical): " << empirical << "\n";
        cout << "  threshold strategy (theory)   : " << theoretical << "\n";
        cout << "  random guess baseline         : " << randomGuess << "\n";
        cout << "  always guess larger baseline  : " << alwaysLarger << "\n\n";
    }

    return 0;
}