using System;

namespace Amphora.RandomComparisons
{
    /// <summary>
    /// Finds the length of the shortest linear feedback shift register that makes a bit sequence, by the
    /// Berlekamp-Massey algorithm over GF(2), with the register and the recent bits packed into words.
    /// </summary>
    public sealed class BerlekampMassey
    {
        private readonly int words;
        private readonly ulong[] connection;
        private readonly ulong[] previous;
        private readonly ulong[] saved;
        private readonly ulong[] recent;

        public BerlekampMassey(int maximumLength)
        {
            words = (maximumLength + 63) / 64 + 1;
            connection = new ulong[words];
            previous = new ulong[words];
            saved = new ulong[words];
            recent = new ulong[words];
        }

        private static int ParityOf(ulong value)
        {
            value ^= value >> 32;
            value ^= value >> 16;
            value ^= value >> 8;
            value ^= value >> 4;
            value ^= value >> 2;
            value ^= value >> 1;
            return (int)(value & 1);
        }

        private void PushRecent(byte bit)
        {
            for (int i = words - 1; i > 0; i--)
            {
                recent[i] = (recent[i] << 1) | (recent[i - 1] >> 63);
            }

            recent[0] = (recent[0] << 1) | bit;
        }

        private void XorPreviousShifted(int shift)
        {
            int wordShift = shift >> 6;
            int bitShift = shift & 63;
            for (int i = words - 1; i >= wordShift; i--)
            {
                ulong value = previous[i - wordShift] << bitShift;
                if (bitShift != 0 && i - wordShift - 1 >= 0)
                {
                    value |= previous[i - wordShift - 1] >> (64 - bitShift);
                }

                connection[i] ^= value;
            }
        }

        /// <summary>
        /// Returns the linear complexity of the first <paramref name="length"/> entries, each 0 or 1.
        /// </summary>
        public int Complexity(byte[] bits, int length)
        {
            Array.Clear(connection, 0, words);
            Array.Clear(previous, 0, words);
            Array.Clear(recent, 0, words);
            connection[0] = 1;
            previous[0] = 1;
            int complexity = 0;
            int lastChange = -1;
            for (int n = 0; n < length; n++)
            {
                PushRecent(bits[n]);
                ulong overlap = 0;
                int limit = complexity / 64 + 1;
                for (int i = 0; i < limit; i++)
                {
                    overlap ^= connection[i] & recent[i];
                }

                if (ParityOf(overlap) == 0)
                {
                    continue;
                }

                Array.Copy(connection, saved, words);
                XorPreviousShifted(n - lastChange);
                if (2 * complexity <= n)
                {
                    complexity = n + 1 - complexity;
                    lastChange = n;
                    Array.Copy(saved, previous, words);
                }
            }

            return complexity;
        }
    }

    /// <summary>
    /// Takes one bit position of consecutive draws in blocks and finds each block's linear complexity.  A random block
    /// of M bits needs a register about M / 2 long, with a known spread.  A generator that is linear in that bit needs
    /// one no longer than its state.
    /// </summary>
    public static class LinearComplexity
    {
        /// <summary>
        /// The chance of each of seven ranges of the complexity's distance from its mean, from the NIST statistical
        /// test suite.
        /// </summary>
        public static readonly double[] Probabilities =
        {
            1.0 / 96, 1.0 / 32, 1.0 / 8, 1.0 / 2, 1.0 / 4, 1.0 / 16, 1.0 / 48,
        };

        /// <summary>
        /// Returns which of the seven ranges a complexity falls in, for a block of an even length.
        /// </summary>
        public static int Category(int complexity, int blockLength)
        {
            double mean = blockLength / 2.0 + 4.0 / 9;
            double distance = complexity - mean + 2.0 / 9;
            if (distance <= -2.5)
            {
                return 0;
            }

            if (distance <= -1.5)
            {
                return 1;
            }

            if (distance <= -0.5)
            {
                return 2;
            }

            if (distance <= 0.5)
            {
                return 3;
            }

            if (distance <= 1.5)
            {
                return 4;
            }

            return distance <= 2.5 ? 5 : 6;
        }

        public static Outcome Run(WordStream stream, Window window, int blockLength, int blocks)
        {
            var machine = new BerlekampMassey(blockLength);
            var bits = new byte[blockLength];
            var observed = new double[Probabilities.Length];
            double totalComplexity = 0;
            for (int block = 0; block < blocks; block++)
            {
                for (int i = 0; i < blockLength; i++)
                {
                    bits[i] = (byte)Windows.EdgeBit(stream.Next(), window);
                }

                int complexity = machine.Complexity(bits, blockLength);
                totalComplexity += complexity;
                observed[Category(complexity, blockLength)]++;
            }

            var expected = new double[Probabilities.Length];
            for (int i = 0; i < expected.Length; i++)
            {
                expected[i] = Probabilities[i] * blocks;
            }

            double degrees = Stats.ChiSquareOfCategories(observed, expected, 5, out double statistic);
            return new Outcome(Stats.ChiSquareTwoSided(statistic, degrees), "bit " + (window == Window.High ? 63 : 0) +
                ", blocks of " + blockLength + ", mean complexity " + (totalComplexity / blocks).ToString("0.0") +
                " against about " + (blockLength / 2.0).ToString("0.0"));
        }
    }
}