using System;

namespace Amphora.RandomComparisons
{
    /// <summary>
    /// Marsaglia's birthday spacings.  Birthdays are drawn on a large calendar and sorted, and the gaps between them
    /// are sorted in turn.  The number of equal gaps follows a Poisson distribution, unless the generator puts its
    /// points on a lattice, which makes equal gaps far more common.
    /// </summary>
    public static class BirthdaySpacings
    {
        public static Outcome Run(WordStream stream, Window window, int dayBits, int birthdays, int samples)
        {
            double mean = Math.Pow(birthdays, 3) / (4.0 * Math.Pow(2, dayBits));
            int categories = (int)(mean * 3 + 20);
            var observed = new double[categories];
            var days = new ulong[birthdays];
            var spacings = new ulong[birthdays];
            ulong calendar = 1UL << dayBits;
            double totalEqualGaps = 0;
            for (int sample = 0; sample < samples; sample++)
            {
                for (int i = 0; i < birthdays; i++)
                {
                    days[i] = Windows.Take(stream.Next(), window, dayBits);
                }

                Array.Sort(days);
                for (int i = 0; i < birthdays - 1; i++)
                {
                    spacings[i] = days[i + 1] - days[i];
                }

                // The calendar is a circle, so the last gap runs around to the first birthday.
                spacings[birthdays - 1] = days[0] + calendar - days[birthdays - 1];
                Array.Sort(spacings);
                int equalGaps = 0;
                for (int i = 1; i < birthdays; i++)
                {
                    if (spacings[i] == spacings[i - 1])
                    {
                        equalGaps++;
                    }
                }

                totalEqualGaps += equalGaps;
                observed[Math.Min(equalGaps, categories - 1)]++;
            }

            double[] expected = Stats.PoissonProbabilities(mean, categories);
            double covered = 0;
            for (int i = 0; i < categories; i++)
            {
                covered += expected[i];
                expected[i] *= samples;
            }

            expected[categories - 1] += (1 - covered) * samples;
            double degrees = Stats.ChiSquareOfCategories(observed, expected, 5, out double statistic);
            return new Outcome(Stats.ChiSquareTwoSided(statistic, degrees), "2^" + dayBits + " days, " + birthdays +
                " birthdays, mean equal gaps " + (totalEqualGaps / samples).ToString("0.000") + " against " +
                mean.ToString("0.###"));
        }
    }
}