using System;

namespace Amphora.RandomComparisons
{
    /// <summary>
    /// Fills square binary matrices from consecutive draws and counts how often each rank comes up.  A generator whose
    /// bits are linear over GF(2) makes rows that depend on each other more or less often than random rows do.
    /// </summary>
    public static class BinaryRank
    {
        /// <summary>
        /// Returns the rank over GF(2) of a matrix of up to 64 rows, each a word whose low bits are the columns.  The
        /// rows are used up.
        /// </summary>
        public static int Rank(ulong[] rows, int size)
        {
            int rank = 0;
            for (int column = 0; column < size; column++)
            {
                ulong mask = 1UL << column;
                int pivot = -1;
                for (int row = rank; row < size; row++)
                {
                    if ((rows[row] & mask) != 0)
                    {
                        pivot = row;
                        break;
                    }
                }

                if (pivot < 0)
                {
                    continue;
                }

                ulong pivotRow = rows[pivot];
                rows[pivot] = rows[rank];
                rows[rank] = pivotRow;
                for (int row = rank + 1; row < size; row++)
                {
                    if ((rows[row] & mask) != 0)
                    {
                        rows[row] ^= pivotRow;
                    }
                }

                rank++;
            }

            return rank;
        }

        /// <summary>
        /// Counts the matrices of a full rank, one short, two short, and further short.  A size of 64 uses whole words.
        /// A size of 32 uses the window's half of each word.
        /// </summary>
        public static Outcome Run(WordStream stream, Window window, int size, long matrices)
        {
            var rows = new ulong[size];
            var observed = new double[4];
            for (long matrix = 0; matrix < matrices; matrix++)
            {
                for (int row = 0; row < size; row++)
                {
                    ulong word = stream.Next();
                    rows[row] = size == 64 ? word : Windows.Take(word, window, size);
                }

                observed[Math.Min(size - Rank(rows, size), 3)]++;
            }

            var expected = new double[4];
            for (int shortBy = 0; shortBy < 3; shortBy++)
            {
                expected[shortBy] = Stats.BinaryRankProbability(size, size - shortBy) * matrices;
            }

            expected[3] = matrices - expected[0] - expected[1] - expected[2];
            double degrees = Stats.ChiSquareOfCategories(observed, expected, 5, out double statistic);
            return new Outcome(Stats.ChiSquareTwoSided(statistic, degrees), size + "x" + size + " matrices of rank " +
                size + ", " + (size - 1) + ", " + (size - 2) + " or less: " + observed[0] + ", " + observed[1] + ", " +
                observed[2] + ", " + observed[3] + " (expected " + expected[0].ToString("0") + ", " +
                expected[1].ToString("0") + ", " + expected[2].ToString("0") + ", " + expected[3].ToString("0") + ")");
        }
    }
}