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 eiher 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.
    */
    public class TorusGridGraphTours
    {

        /// <summary>
        /// Computes the number of directed Hamilton paths on the (n,n)-torus grid graph. (A137891)
        /// </summary>
        /// <remarks>
        /// The (n,m)-torus grid graph is the disjoint graph join of Cn and Cm. (has n+m vertices).
        /// </remarks>
        public static BigInteger TorusGridGraphHamiltonPathCount( int n )
        {
            return TorusGridGraphHamiltonPathCount( n, n );
        }

        /// <summary>
        /// Computes the number of directed Hamilton paths on the (n,m)-torus grid graph.
        /// </summary>
        /// <remarks>
        /// The (n,m)-torus grid graph is the disjoint graph join of Cn and Cm. (has n+m vertices).
        /// </remarks>
        public static BigInteger TorusGridGraphHamiltonPathCount( int n, int m )
        {

            // The essence of the computation is to loop over all pairs of partitions that have the same number of pieces
            // or where the first partition has one more piece than the second. (Each path through the graph
            // corresponds to a zig-zag from one side to the other where the vertices are separated into groups of n and m
            // in the most obvious manner. The number of hops must balance if the path does not
            // return to the original side, otherwise there will be an extra hop back to the original side).
            // For each pair of partitions we form a product of all the possibilities and sum to get the total.
            // The algorithm below does the computation in two phases to avoid order n squared processing,
            // which improves efficiency when there are many partitions. This also leads a long way towards a formula.

            BigInteger[] factorials = Factorials().Take( Math.Min( m, n ) + 2 ).ToArray();

            Func<IReadOnlyList<BigInteger>, IReadOnlyList<BigInteger>, BigInteger> CrossMultiply = ( v1, v2 ) =>
            {
                var qry = from partCount in Enumerable.Range( 1, Math.Min( v1.Count - 1, v2.Count ) )
                          select v1[ partCount ] * factorials[ partCount ] * factorials[ partCount - 1 ]
                                 * ((partCount < v2.Count ? partCount * v2[ partCount ] : 0) + v2[ partCount - 1 ]);
                return Sum( qry );
            };
            var pn = ComputePartitionVector( n );
            var pm = m == n ? pn : ComputePartitionVector( m );
            return CrossMultiply( pn, pm ) + CrossMultiply( pm, pn );

        }

        /// <summary>
        /// Computes the number of Hamilton cycles on the (n,n)-torus grid graph.
        /// </summary>
        public static BigInteger TorusGridGraphHamiltonCycleCount( int n )
        {
            return TorusGridGraphHamiltonCycleCount( n, n );
        }

        /// <summary>
        /// Computes the number of Hamilton cycles on the (n,m)-torus grid graph.
        /// </summary>
        public static BigInteger TorusGridGraphHamiltonCycleCount( int n, int m )
        {
            BigInteger[] factorials = Factorials().Take( Math.Min( m, n ) + 1 ).ToArray();

            Func<IReadOnlyList<BigInteger>, IReadOnlyList<BigInteger>, BigInteger> CrossMultiply = ( v1, v2 ) =>
            {
                var qry = from partCount in Enumerable.Range( 1, Math.Min( v1.Count, v2.Count ) - 1 )
                          select v1[ partCount ] * v2[ partCount ] * factorials[ partCount ] * factorials[ partCount - 1 ];
                return Sum( qry );
            };

            var pn = ComputePartitionVector( n );
            var pm = m == n ? pn : ComputePartitionVector( m );
            return CrossMultiply( pn, pm ) / 2;
        }

        /// <summary>
        /// Computes intermediate summary information on all partitions of n broken out by number of parts.
        /// </summary>
        private static IReadOnlyList<BigInteger> ComputePartitionVector( int n )
        {
            // 3 different implementations to choose from. (uncomment one)
            return ComputePartitionVectorCombinatorial( n );
            //return ComputePartitionVectorAlternate( n );
            //return ComputePartitionVectorFormula( n );
        }

        /// <summary>
        /// Combinatorial implementation.
        /// </summary>
        private static IReadOnlyList<BigInteger> ComputePartitionVectorCombinatorial( int n )
        {
            BigInteger[] result = new BigInteger[ n + 1 ];

            foreach (var partition in Partition.Enumerate( n, n, p => p ))
            {
                result[ PartCount( partition ) ] += FullArrangementFactor( partition );
            }
            return result;
        }

        /// <summary>
        /// Equivalent to the above. (probably slower)
        /// </summary>
        private static IReadOnlyList<BigInteger> ComputePartitionVectorAlternate( int n )
        {
            BigInteger[] result = new BigInteger[ n + 1 ];
            result[ 1 ] = n <= 2 ? n : 2 * n;
            for (int initial = 1; initial < n; ++initial)
            {
                foreach (var partition in Partition.Enumerate( n - initial, n, p => p ))
                {
                    result[ PartCount( partition ) + 1 ] += SimpleArrangementFactor( partition ) * (initial != 1 ? 2 * initial : 1);
                }
            }
            return result;
        }

        /// <summary>
        /// Equivalent to the above.
        /// </summary>
        private static IReadOnlyList<BigInteger> ComputePartitionVectorFormula( int n )
        {
            if (n == 2)
            {
                // handle this as a special case to avoid complicating formula
                return new BigInteger[] { 0, 2, 1 };
            }

            BigInteger[] factorials = Factorials().Take( n + 1 ).ToArray();

            Func<int, int> A130130 = t => t >= 2 ? 2 : t;

            Func<int, int, BigInteger> A266213 = ( t, r ) =>
            {
                // Sum_ { k=0..min(n,r)} binomial( r-1, k-1)*binomial( n, k)*2^k.
                var qry = from k in Enumerable.Range( 1, Math.Max( 0, Math.Min( t, r ) ) )
                          select factorials[ r - 1 ] * factorials[ t ] * BigInteger.Pow( 2, k )
                            / (factorials[ k - 1 ] * factorials[ r - k ] * factorials[ k ] * factorials[ t - k ]);
                return r == 0 ? 1 : Sum( qry );
            };


            //Sum_ { k=1..n } 2*k!*b(n,k)*(k!*b(n,k)+(k-1)!*b(n,k-1)) 
            // where b(n,0)=0, b(n,k)=Sum_{ j=1..n-k+1 } j*A130130(j)*A266213(k-1,n-j-k+1) for k>0, n!=2.
            BigInteger[] result = new BigInteger[ n + 1 ];
            for (int k = 1; k <= n; ++k)
            {
                for (int initial = 1; initial <= n - k + 1; ++initial)
                {
                    result[ k ] += initial * A130130( initial ) * A266213( k - 1, n - initial - k + 1 );
                }
            }

            return result;
        }

        /// <summary>
        /// Returns number of parts in a partition.
        /// </summary>
        private static int PartCount( IReadOnlyList<Tuple<int, int>> partitionPieces )
        {
            return partitionPieces.Sum( piece => piece.Item2 );
        }

        /// <summary>
        /// Computes required multiplication factor for a given partition.
        /// </summary>
        private static BigInteger FullArrangementFactor( IReadOnlyList<Tuple<int, int>> pieces )
        {
            int partCount = PartCount( pieces );
            int unitCount = (short)(pieces[ pieces.Count - 1 ].Item1 == 1 ? pieces[ pieces.Count - 1 ].Item2 : 0);

            if (partCount > 1 || pieces[ pieces.Count - 1 ].Item1 != 2)
            {
                BigInteger[] factorials = Factorials().Take( partCount + 1 ).ToArray();
                BigInteger numerator = 0;
                BigInteger denominator = 1;
                foreach (var piece in pieces)
                {
                    numerator += piece.Item1 * piece.Item2;
                    denominator *= factorials[ piece.Item2 ];
                }
                numerator *= factorials[ partCount - 1 ];
                if (numerator % denominator != 0)
                {
                    throw new ApplicationException( "Internal error - division should be exact." );
                }
                return numerator / denominator * BigInteger.Pow( 2, partCount - unitCount );
            }
            else
            {
                // this is a special case that uniquely applies to partition of 2 into 1 part of 2.
                return 2;
            }
        }

        /// <summary>
        /// A 'reduced' version of the above that does not treat the initial piece as a special case.
        /// </summary>
        private static BigInteger SimpleArrangementFactor( IReadOnlyList<Tuple<int, int>> pieces )
        {

            int pieceCount = (short)pieces.Sum( piece => piece.Item2 );
            int unitCount = (short)(pieces[ pieces.Count - 1 ].Item1 == 1 ? pieces[ pieces.Count - 1 ].Item2 : 0);
            BigInteger[] factorials = Factorials().Take( pieceCount + 1 ).ToArray();
            // compute multinomial of pieces
            BigInteger denominator = 1;
            foreach (var piece in pieces)
            {
                denominator *= factorials[ piece.Item2 ];
            }
            return factorials[ pieceCount ] / denominator * BigInteger.Pow( 2, pieceCount - unitCount );
        }


        /// <summary>
        /// Helper to sum up a sequence of integers.
        /// </summary>
        private static BigInteger Sum( IEnumerable<BigInteger> values )
        {
            BigInteger total = 0;
            foreach (var v in values)
            {
                total += v;
            }
            return total;
        }

        /// <summary>
        /// Helper to get factorials. Returns zero based sequence.
        /// </summary>
        private static IEnumerable<BigInteger> Factorials()
        {
            BigInteger product = 1;
            for (int m = 1; ; ++m)
            {
                yield return product;
                product *= m;
            }
        }

    }

    /// <summary>
    /// Class for working with partitions.
    /// </summary>
    public static class Partition
    {
        /// <summary>
        /// Gets the number of partitions of a given integer into parts of no more than a given maximum size.
        /// </summary>
        public static int PartitionCount( int total, int maxSize )
        {
            return Enumerate( total, maxSize, p => 1 ).Count();
        }

        /// <summary>
        /// Enumerates all partitions of a given integer into parts of no more than a given maximum size.
        /// </summary>
        /// <remarks>
        /// Each generated partition consits of parts ordered in descending order of size and is passed to the capture
        /// function for processing. Item1 in each tuple is the partition size and Item2 is the number of pieces of that size.
        /// </remarks>
        public static IEnumerable<T> Enumerate<T>( int total, int max, Func<IReadOnlyList<Tuple<int, int>>, T> capture )
        {
            if (total <= 0)
            {
                throw new ArgumentException( "total", "Must be a positive integer." );
            }
            if (max <= 0)
            {
                throw new ArgumentException( "max", "Must be a positive integer." );
            }

            if (total > 0)
            {
                if (max > total)
                {
                    max = total;
                }
                var values = new List<Tuple<int, int>>();
                do
                {
                    if (total > 0)
                    {
                        int d = total / max;
                        values.Add( Tuple.Create( max, d ) );
                        total -= max * d;
                        max = total;
                    }
                    if (total == 0)
                    {
                        yield return capture( values );
                        int k = values.Count;
                        var top = values[ --k ];
                        if (top.Item1 == 1)
                        {
                            total += top.Item2;
                            values.RemoveAt( k );
                            --k;
                        }
                        if (k >= 0)
                        {
                            top = values[ k ];
                            total += top.Item1;
                            max = top.Item1 - 1;
                            if (top.Item2 > 1)
                            {
                                values[ k ] = Tuple.Create( top.Item1, top.Item2 - 1 );
                            }
                            else
                            {
                                values.RemoveAt( k );
                            }
                        }
                    }

                } while (max != 0);
            }
        }
    }

}