using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;



namespace Application
{

    /* Author: Andrew Howroyd. (c) 2016
       Terms of use: You may freely use this software for personal, open source, educational or commercial use.
       You may release this software under either BSD or GPLL open source licences.
       If you modify, adapt or reformat this code, please state clearly that the code has been modified from the original.
       No waranty.
    */

    /// <summary>
    /// Class for computing the number of matchings or perfect matchings in the rook's graph K_n x K_n.
    /// It can also be used to compute the number of matchings in a bishop's graph or other induced subgraph
    /// defined by a Ferrers diagram.
    /// The algorithm is straight forward, and illustrates many of the techniques needed to compute
    /// the number of Hamilton cycles or paths in these graphs. (details more complicated, but conceptually similar).
    /// </summary>
    public class LatticeMatchings
    {

        /// <summary>
        /// Computes the number of perfect matchings in the m X n lattice graph.
        /// </summary>
        public static BigInteger LatticePerfectMatchingCount( int width, int height )
        {
            var processor = new LatticeMatchings( RectangularHeights( width, height ), true );
            List<KeyValuePair<BigInteger, BigInteger>> counts = processor.InitialCounts().ToList();
            return ProcessLayers( processor, counts );
        }

        /// <summary>
        /// Computes the number of matchings in the m x n lattice graph.
        /// </summary>
        public static BigInteger LatticeMatchingCount( int width, int height )
        {
            var processor = new LatticeMatchings( RectangularHeights( width, height ), false );
            List<KeyValuePair<BigInteger, BigInteger>> counts = processor.InitialCounts().ToList();
            return ProcessLayers( processor, counts );
        }

        /// <summary>
        /// Heights for a rectangular array.
        /// </summary>
        public static IEnumerable<int> RectangularHeights( int width, int height )
        {
            return Enumerable.Repeat( height, width );
        }

        /// <summary>
        /// Top level function for computing the total number of matchings.
        /// </summary>
        private static BigInteger ProcessLayers( LatticeMatchings processor, List<KeyValuePair<BigInteger, BigInteger>> counts )
        {
            while (!processor.IsFinal)
            {
                processor = processor.GetNextLayer();
                counts = processor.Process( counts ).ToList();
                counts.TrimExcess();
            }
            processor = null;
            BigInteger total = BigInteger.Zero;
            foreach (var kv in counts)
            {
                total += (BigInteger)kv.Value;
            }
            return total;
        }

        private readonly int[] m_Heights;
        private readonly int m_Width;
        private readonly bool m_Perfect;
        private readonly int m_StartColumn;
        private readonly int[] m_LinkCounts;
        private int m_CrossLinks;
        private Dictionary<BigInteger, BigInteger> m_States;
        private long m_Transitions;

        /// <summary>
        /// Constructor.
        /// </summary>
        public LatticeMatchings( IEnumerable<int> heights, bool perfect ) : this( heights, perfect, 0 )
        {

        }

        /// <summary>
        /// Internal constructor.
        /// </summary>
        private LatticeMatchings( IEnumerable<int> heights, bool perfect, int startColumn )
        {
            m_Heights = heights.Where( h => h >= 0 ).OrderBy( h => h ).ToArray();
            int width = m_Heights.Length;
            m_Width = width;
            m_LinkCounts = new int[ width + 1 ];
            m_StartColumn = startColumn < 0 ? width - 1 : startColumn;
            m_Perfect = perfect;
        }

        /// <summary>
        /// If this is the final layer.
        /// </summary>
        public bool IsFinal
        {
            get { return m_Heights[ m_Width - 1 ] == 0 && m_StartColumn == 0; }
        }

        /// <summary>
        /// Gets the next processing layer.
        /// </summary>
        public LatticeMatchings GetNextLayer()
        {
            return new LatticeMatchings( m_StartColumn > 0 ? m_Heights : m_Heights.Select( h => h - 1 ), m_Perfect, m_StartColumn - 1 );
        }

        /// <summary>
        /// Initial state/count list.
        /// </summary>
        public IEnumerable<KeyValuePair<BigInteger, BigInteger>> InitialCounts()
        {
            var states = new Dictionary<BigInteger, BigInteger>();
            m_States = states;
            Accumulate( PackState(), BigInteger.One );
            m_States = null;
            return states;
        }

        /// <summary>
        /// Processes a single cell of the lattice. Input and output are lists of state/count pairs.
        /// </summary>
        public IEnumerable<KeyValuePair<BigInteger, BigInteger>> Process( IEnumerable<KeyValuePair<BigInteger, BigInteger>> input )
        {
            var states = new Dictionary<BigInteger, BigInteger>();
            m_States = states;
            foreach (var kv in input)
            {
                UnpackState( kv.Key );
                ExploreOptions( m_StartColumn, kv.Value );
            }
            m_States = null;
            return states;
        }

        private void ExploreOptions( int col, BigInteger multiplicity )
        {
            int n = m_LinkCounts[ col ];
            if (n > 0)
            {
                // close a vertica link
                m_LinkCounts[ col ] = n - 1;
                MergeResults( col, multiplicity * n );
                m_LinkCounts[ col ] = n;
            }
            if (n < m_Heights[ col ])
            {
                // open a vertical link
                m_LinkCounts[ col ] = n + 1;
                MergeResults( col, multiplicity );
                m_LinkCounts[ col ] = n;
            }
            if (n <= m_Heights[ col ])
            {
                if (m_CrossLinks > 0)
                {
                    // close a horizontal link
                    --m_CrossLinks;
                    MergeResults( col, multiplicity * (m_CrossLinks + 1) );
                    ++m_CrossLinks;
                }
                if (m_CrossLinks < col)
                {
                    // open a horizontal link
                    ++m_CrossLinks;
                    MergeResults( col, multiplicity );
                    --m_CrossLinks;
                }
                if (!m_Perfect)
                {
                    // leave vertex unconnected
                    MergeResults( col, multiplicity );
                }
            }
        }

        private void MergeResults( int col, BigInteger multiplicity )
        {
            // prune states that cannot be completed due to competition between vertical and horizontal links.
            int t = m_CrossLinks;
            for (int p = 0; p < col && t > 0; ++p)
            {
                if (m_LinkCounts[ p ] <= m_Heights[ p ])
                {
                    --t;
                }
            }
            if (t == 0)
            {
                // maintain sort order of link counts. This reduces the number of states
                int h = m_Heights[ col ];
                int v = m_LinkCounts[ col ];
                int p = col + 1;

                while (p < m_Width && m_Heights[ p ] == h && m_LinkCounts[ p ] > v)
                {
                    m_LinkCounts[ p - 1 ] = m_LinkCounts[ p ];
                    ++p;
                }
                m_LinkCounts[ p - 1 ] = v;

                // add to output dictionary.
                Accumulate( PackState(), multiplicity );

                // need to restore order to orignal
                while (--p > col)
                {
                    m_LinkCounts[ p ] = m_LinkCounts[ p - 1 ];
                }
                m_LinkCounts[ col ] = v;
            }
        }

        /// <summary>
        /// Returns the number of transitions processed. Gives a measure of work done.
        /// </summary>
        public long TransitionCount
        {
            get { return m_Transitions; }
        }

        /// <summary>
        /// Add state and count to current output dictionary.
        /// </summary>
        private void Accumulate( BigInteger state, BigInteger multiplicity )
        {
            if (!multiplicity.IsZero)
            {
                BigInteger n;
                m_States[ state ] = m_States.TryGetValue( state, out n ) ? n + multiplicity : multiplicity;
                ++m_Transitions;
            }
        }

        /// <summary>
        /// Saves current state into a BigInteger word.
        /// </summary>
        private BigInteger PackState()
        {
            BigInteger state = 0;
            m_LinkCounts[ m_Width ] = m_StartColumn > 0 ? m_CrossLinks : 0;
            for (int col = 0; col < m_LinkCounts.Length; ++col)
            {
                int v = m_LinkCounts[ col ];
                while (v >= 3)
                {
                    state = (state << 2) | 3;
                    v -= 3;
                }
                state = (state << 2) | v;
            }
            return state;
        }

        /// <summary>
        /// Restores current state.
        /// </summary>
        private void UnpackState( BigInteger state )
        {

            for (int col = m_LinkCounts.Length; --col >= 0;)
            {
                int v = (int)(state & 3);
                state >>= 2;
                while ((state & 3) == 3)
                {
                    v += 3;
                    state >>= 2;
                }
                m_LinkCounts[ col ] = v;
            }
            m_CrossLinks = m_LinkCounts[ m_Width ];
        }
    }
}