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

namespace Application
{

    /* Author: Andrew Howroyd. (c) 2015
       Terms of use: You may freely use this software for personal, open source, educational or commercial use.
       However, you may not redistribute or use this software in conjunction with any modification or improvement for which use is
       not fully granted to all under the terms of this or GPL/LGPL license. (to use freely you must be willing to give freely)
       You may release this code under the GNU Lesser Public License (LGPL v2.1: http://www.gnu.org/licenses/old-licenses/lgpl-2.1.html)
       If you modify, adapt or reformat this code, please state clearly that the code has been modified from the original.
       No waranty.
    */

    /*
        This software is not a complete program. It is left to the end user to provide application code, api and top level enumeration functions.
        It is more intended as a starter tool-kit for those wanting to explore meanders or reproduce some of the known results.
        The code was compiled using Visual Studio Community 2015 and was run on the Windows platform using .Net 4.5.
        For BigInteger, the System.Numerics package must be added as a reference to the project. This is part of the standard .Net library.

        The basic meander transfer enumeration method is implemented in the MeanderProblem class. Sub classes of this provide
        customizations for different sequences. Included here are BasicMeanderProblem, MeandersByComponents and ParallelRoadsMeanderProblem.
        The first two enumerate the most important core sequences, whilst the ParallelRoadsMeanderProblem is a demonstration
        of how the same algorithm can be extended to a more complex problem - by pretending the multiple roads are a single
        snake like road and adding book keeping to track the meander from one bend to the next.
        The algorithm has also been successfully used to enumerate a variety of other meander sequences. (A060148, A077056, A060972, A259974).
        Code for these is not included here as it would add unnecessary clutter with little additional insight.

        In addition to a sub-class of BasicMeanderProblem a state machine processor is required. Two are included here.
        SimpleProcessor is the most natutal and simple implementation, but will eventually fail at about 40 bridges due to
        insufficient memory to contain the entire state space. This is suitable for many needs.
        DivideAndConquerProcessor is a processor that reduces memory use at the cost of additional cpu by sub-dividing
        the problem into multiple partitions which are processed independently.

        The method is generally applicable to any meander problem that can be reduced to a simple left to right scan of the meander.
        More generally the 'transfer-matrix' method is applicable to a very wide class of enumeration problem.

        Sample usage:
        
        public BigInteger OpenMeanderCount( int length )
        {
            int nn = length - 1;
            var processor = new SimpleProcessor<BigInteger>() { CreateStateMachine = k => new BasicMeanderProblem( k ) };
            return processor.Process( nn, new BasicMeanderProblem( nn ).OpenMeanderInitialStates );
        }

        ArchConfiguration which is un-related and indpendent of the other classes is a set of utilities for working
        with arch configuration and menander permutations.   
        
    */

    #region MeanderProblem
    /// <summary>
    /// Base class for meander enumeration problems.
    /// </summary>
    /// <remarks>
    /// Algorithm: 
    /// Based on: I. Jensen, A transfer matrix approach to the enumeration of plane meanders, J. Phys. A 33 (2000), no. 34, 5953-5963. 
    /// http://arXiv.org/abs/cond-mat/0008178
    /// There are a few minor implementation differences.
    ///   - Two separate bit sequences are used for the lower and upper segments, with the least significant bits closest to the road.
    ///     These are inter-woven into a single integer valued state.
    ///   - Bit value 1 is used to indicate that an arch originates on this side of the road and 0 indicates on the other side. This gives better symmetry when processing.
    ///   - This algorithm does not have any notion of a free end - there is always a loop. Open meanders can be enumerated with an alternative initial state.
    /// 
    /// Performance:
    /// The primary limiting factor of the algorithm as presented here is memory rather than cpu. Approximately 2GB is required to get to about 40 bridges.
    /// However, by adopting a multi-stage process it is possible to some extent to trade memory usage for additional cpu.
    /// The algorithm as presented here uses BigInteger for both states and counts. For better cpu performance,  64-bit long integers can be used
    /// at least for the open meander problem up to about 60 bridges before overflow will be an issue. 
    /// </remarks>
    abstract class MeanderProblem
    {
        protected const int WORD_SHIFT = 2;
        protected const int WORD_MASK = (1 << WORD_SHIFT) - 1;
        protected readonly int m_LayerIndex;

        protected enum Action
        {
            None, NewArch, JoinArch, MoveUp, MoveDown, CloseLoop
        }

        /// <summary>
        /// Constructor.
        /// </summary>
        protected MeanderProblem( int layerIndex )
        {
            m_LayerIndex = layerIndex;
        }

        /// <summary>
        /// Index of the transition layer. Useful for diagnostics/progress reporting.
        /// </summary>
        public int LayerIndex
        {
            get { return m_LayerIndex; }
        }

        /// <summary>
        /// Enumerate next states from previous state.
        /// </summary>
        protected IEnumerable<T> EnumeratePossibilities<T>( BigInteger state, Func<Action, BigInteger, BigInteger, T> capture )
        {
            // split the incoming state into lower and upper halves for ease of processing.
            var lower = ExtractLower( state );
            var upper = (state - lower) >> 1;

            // There are four possible state transitions.
            // 1. Leg of arch crosses from below road to above road. Allowed if not at edge guard.
            if (!lower.IsOne)
            {
                yield return capture( Action.MoveUp, lower >> WORD_SHIFT, (upper << WORD_SHIFT) ^ CrossRoad( lower ) );
            }
            // 2. Leg of arch crosses from above road to below road. Allowed if not at edge guard.
            if (!upper.IsOne)
            {
                yield return capture( Action.MoveDown, (lower << WORD_SHIFT) ^ CrossRoad( upper ), upper >> WORD_SHIFT );
            }
            // 3. Introduction of new arch. Allowed always.
            yield return capture( Action.NewArch, (lower << WORD_SHIFT) | BigInteger.One, (upper << WORD_SHIFT) | BigInteger.One );

            // 4. Arch connection. Possible if not at edge guard on either side.
            if (!lower.IsOne && !upper.IsOne)
            {
                var test = ClosingAction( lower, upper );
                if (test != Action.None)
                {
                    yield return capture( test, JoinArch( lower, upper ) >> WORD_SHIFT, JoinArch( upper, lower ) >> WORD_SHIFT );
                }
            }
        }

        protected virtual Action ClosingAction( BigInteger lower, BigInteger upper )
        {
            return lower.IsEven || upper.IsEven ? Action.JoinArch : Action.CloseLoop;
        }

        /// <summary>
        /// When crossing the road, we need to invert the bit that indicates if the arch originated on this side of the road.
        /// </summary>
        protected static BigInteger CrossRoad( BigInteger v )
        {
            return v.IsEven ? BigInteger.One : BigInteger.Zero;
        }

        /// <summary>
        ///  Handles the mechanics of re-joining arches.
        /// </summary>
        private static BigInteger JoinArch( BigInteger v, BigInteger alt )
        {
            // see Iwan Jenson paper to understand what is going on here. (with illustrations)
            if (v.IsEven && !alt.IsEven)
            {
                int n = 0;
                BigInteger bit = BigInteger.One;
                while (n >= 0)
                {
                    bit <<= WORD_SHIFT;
                    n += (v & bit).IsZero ? 1 : -1;
                }
                v ^= bit;
            }
            return v;
        }

        /// <summary>
        /// Test for enclosure. Used to identitfy forest meander systems. (A060148)
        /// </summary>
        /// <param name="v">Upper or lower half (doesn't matter which)</param>
        protected bool IsEnclosed( BigInteger v )
        {
            bool enclosed = true;
            for (var bit = (BigInteger)2; bit <= v; bit <<= WORD_SHIFT)
            {
                enclosed = !enclosed;
            }
            return enclosed;
        }

        /// <summary>
        /// Combines the states for the lower and upper halves into a single integer.
        /// Each half only uses even numbered bits, so they can be interleaved with a simple bit shift.
        /// </summary>
        protected static BigInteger Pack( BigInteger lower, BigInteger upper )
        {
            if (lower.IsZero || upper.IsZero)
            {
                throw new ApplicationException( "invalid state" );
            }
            return lower | (upper << 1);
        }

        /// <summary>
        /// Variation on Pack that may swap lower and upper words. This reduces the total
        /// number of states by about half, but is not safe in those cases where the distinction
        /// between upper and lower is important. (This is purely a questionable optimisation
        /// that in a few cases allows an extra term to be computed)
        /// </summary>
        protected static BigInteger PackSymmetrical( BigInteger lower, BigInteger upper )
        {
            return lower < upper ? (upper | (lower << 1)) : (lower | (upper << 1));
        }

        /// <summary>
        /// Extracts the lower arch configuration from an encoding.
        /// The upper arch configuration can the be obtained by subtraction.
        /// This undoes the encoding of Pack.
        /// </summary>
        protected static BigInteger ExtractLower( BigInteger v )
        {
            BigInteger mask = 0x5555555555555555;
            int bits = 64;
            while (mask < v)
            {
                mask |= (mask << bits);
                bits += bits;
            }
            return v & mask;
        }


        /// <summary>
        /// The initial state of the traversal.
        /// </summary>
        protected static BigInteger DefaultInitialState
        {
            // The initial state consists of two edge guard bits.
            get { return Pack( BigInteger.One, BigInteger.One ); }
        }

        /// <summary>
        /// Provides a heuristic to group states into partitions. Used by the DivideAndConquerProcessor.
        /// The aim is to put similar states into the same partition, but there is no guarantee that this achieves that goal.
        /// </summary>
        public static int Cluster( BigInteger state, int level )
        {
            Func<BigInteger, BigInteger> Reduce = v =>
            {
                var limit = BigInteger.One << (2 + (WORD_SHIFT * level));
                if (v < limit)
                {
                    return BigInteger.One;
                }
                while (v >= limit)
                {
                    v >>= WORD_SHIFT;
                }
                return v;
            };
            var lower = ExtractLower( state );
            var upper = (state - lower) >> 1;
            return (int)Pack( Reduce( lower ), Reduce( upper ) );
        }

    }
    #endregion

    #region IStateMachine
    /// <summary>
    /// Interface to abstract state enumeration process.
    /// </summary>
    /// <typeparam name="T">Type of state</typeparam>
    interface IStateMachine<T>
    {
        int LayerIndex
        { get; }
        IEnumerable<T> Enumerate( T state );
    }
    #endregion

    #region BasicMeanderProblem
    /// <summary>
    /// Processing component to determine number of open meanders. (A005316).
    /// </summary>
    class BasicMeanderProblem : MeanderProblem, IStateMachine<BigInteger>
    {
        private readonly BigInteger m_Limit;
        private readonly bool m_IsOdd;

        public BasicMeanderProblem( int remainingBridges ) : base( remainingBridges )
        {
            m_Limit = BigInteger.One << (2 + (WORD_SHIFT * (remainingBridges + 1)));
            m_IsOdd = (remainingBridges & 1) == 1;
        }

        /// <summary>
        /// Initial states to enumerate open meansers.
        /// </summary>
        public IEnumerable<BigInteger> OpenMeanderInitialStates
        {
            get
            {
                // by modifying the intial state we can count open meanders using a closed meander algorithm.
                yield return m_IsOdd ? Pack( ((1 << WORD_SHIFT) | 1) << WORD_SHIFT, 1 ) : Pack( (1 << WORD_SHIFT) | 1, (1 << WORD_SHIFT) | 1 );
            }
        }

        /// <summary>
        /// Initial states used to enumerate semi-meanders. (A000682)
        /// </summary>
        public IEnumerable<BigInteger> SemiMeanderInitialStates
        {
            get
            {
                BigInteger bits = m_IsOdd ? 1 : (1 << WORD_SHIFT) | 1;
                var state = Pack( bits, bits );
                while (state < m_Limit)
                {
                    yield return state;
                    bits = (((bits << WORD_SHIFT) | BigInteger.One) << WORD_SHIFT) | BigInteger.One;
                    state = Pack( bits, bits );
                }
            }
        }

        /// <summary>
        /// Initial states for rotationally symmetric loops on two parallel roads. (part of A086031)
        /// gives 1,3,16,105,786,6398,55280,499293,4667290,44840730
        /// </summary>
        public IEnumerable<BigInteger> ParallelRoadSymmetricMeanderInitialStates
        {
            get
            {
                BigInteger bits = (1 << WORD_SHIFT) | 1;
                var state = Pack( bits << WORD_SHIFT, BigInteger.One );
                int r = 1;
                while (state < m_Limit)
                {
                    yield return state;
                    r += 2;
                    bits = (((bits << WORD_SHIFT) | BigInteger.One) << WORD_SHIFT) | BigInteger.One;
                    state = Pack( bits << (WORD_SHIFT * r), BigInteger.One );
                }
            }
        }
        /// <summary>
        /// Enumerate next states from previous state.
        /// </summary>
        public IEnumerable<BigInteger> Enumerate( BigInteger state )
        {
            return EnumeratePossibilities( state, ( action, lower, upper ) => Pack( lower, upper ) )
                     .Where( next => m_Limit.Sign < 0 || next < m_Limit );
        }

        protected override Action ClosingAction( BigInteger lower, BigInteger upper )
        {
            return base.ClosingAction( lower, upper ) == Action.JoinArch ? Action.JoinArch : Action.None;
        }
    }
    #endregion

    #region MeandersByComponents
    /// <summary>
    /// Processing component to count meander systems or semi-menader systems by number of component.
    /// A006657, A008828, A046721, A046726
    /// </summary>
    class MeandersByComponents : MeanderProblem, IStateMachine<Tuple<int, BigInteger>>
    {
        private readonly BigInteger m_Limit;
        private readonly int? m_MaxComponents;
        private readonly int m_RemainingBridges;

        public MeandersByComponents( int remainingBridges, int? maxComponents = null ) : base( remainingBridges )
        {
            m_RemainingBridges = remainingBridges;
            m_Limit = BigInteger.One << (2 + (WORD_SHIFT * remainingBridges));
            m_MaxComponents = remainingBridges == 0 || !maxComponents.HasValue ? maxComponents : maxComponents.Value - 1;
        }

        /// <summary>
        /// Initial states for closed meander systems.
        /// </summary>
        public IEnumerable<Tuple<int, BigInteger>> ClosedMeanderInitialStates
        {
            get { yield return Tuple.Create( 0, DefaultInitialState ); }
        }

        /// <summary>
        /// Initial states for semi-meander sysytems.
        /// </summary>
        public IEnumerable<Tuple<int, BigInteger>> SemiMeanderInitialStates( Func<int, bool> windingPredicate )
        {
            var bits = BigInteger.One;
            var state = Pack( bits, bits );
            int winding = 0;
            while (state < m_Limit)
            {
                if ((winding & 1) == (m_RemainingBridges & 1) && (windingPredicate == null || windingPredicate( winding )))
                {
                    yield return Tuple.Create( 0, state );
                }
                ++winding;
                bits = (bits << WORD_SHIFT) | BigInteger.One;
                state = Pack( bits, bits );
            }
        }

        public IEnumerable<Tuple<int, BigInteger>> Enumerate( Tuple<int, BigInteger> state )
        {
            int n = state.Item1;
            return EnumeratePossibilities( state.Item2,
                  ( action, lower, upper ) => Tuple.Create( action == Action.CloseLoop ? n + 1 : n, PackSymmetrical( lower, upper ) ) )
                .Where( next => (m_Limit.Sign < 0 || next.Item2 < m_Limit) && (!m_MaxComponents.HasValue || state.Item1 <= m_MaxComponents.Value) );
        }
    }
    #endregion

    #region ParallelRoadsMeanderProblem
    /// <summary>
    /// Processing component for meanders with multiple parallel roads. (A076876, A206432, A204352 etc). Can also handle loops on multiple parallel roads (A086031).
    /// </summary>
    class ParallelRoadsMeanderProblem : MeanderProblem, IStateMachine<Tuple<int, BigInteger>>
    {
        private readonly BigInteger m_Limit;
        private readonly int m_MaxRoads;
        private readonly bool m_IsPath;
        private readonly int m_Remaining;
        private readonly int m_Sign;

        public ParallelRoadsMeanderProblem( int remaining, int maxRoads, bool isPath, bool isOdd ) : base( remaining )
        {
            m_Limit = BigInteger.One << (2 + (WORD_SHIFT * (remaining + 1)));
            m_Remaining = remaining;
            m_MaxRoads = maxRoads;
            m_IsPath = isPath;
            m_Sign = isOdd ? 0 : 1;
        }

        public Tuple<int, BigInteger> InitialState
        {
            get
            {
                return Tuple.Create( m_IsPath ? -1 : 1,
              m_Remaining % 2 == 0 ? Pack( (BigInteger.One << WORD_SHIFT) | BigInteger.One, (BigInteger.One << WORD_SHIFT) | BigInteger.One )
              : Pack( BigInteger.One, ((BigInteger.One << WORD_SHIFT) | BigInteger.One) << WORD_SHIFT ) );
            }
        }

        public IEnumerable<Tuple<int, BigInteger>> Enumerate( Tuple<int, BigInteger> state )
        {
            return EnumeratePossibilities( state.Item2,
                    ( action, lower, upper ) => Tuple.Create( state.Item1 < 0 && upper.IsOne ? -state.Item1 : state.Item1, Pack( lower, upper ) ) )
                .Where( next => m_Limit.Sign < 0 || next.Item2 < m_Limit )
                .SelectMany( Fork );
        }

        private IEnumerable<Tuple<int, BigInteger>> Fork( Tuple<int, BigInteger> state )
        {
            // There is always the option to stay on the current road.
            yield return state;

            // If we are not yet at the limit, also look into the possibility of moving to the next.
            if (state.Item1 < m_MaxRoads && -state.Item1 < m_MaxRoads)
            {
                var lower = ExtractLower( state.Item2 );
                var upper = (state.Item2 - lower) >> 1;
                bool canAdvance = (state.Item1 & 1) == m_Sign ? lower.IsOne
                                : state.Item1 < 0 ? upper == ((1 << WORD_SHIFT) | 1) : upper.IsOne;
                if (canAdvance)
                {
                    yield return Tuple.Create( state.Item1 + (state.Item1 < 0 ? -1 : 1), state.Item2 );
                }
            }
        }

        protected override Action ClosingAction( BigInteger lower, BigInteger upper )
        {
            return base.ClosingAction( lower, upper ) == Action.JoinArch ? Action.JoinArch : Action.None;
        }
    }
    #endregion

    #region SimpleProcessor
    /// <summary>
    /// Simple processing engine that will work with any kind of state type.
    /// </summary>
    class SimpleProcessor<StateType>
    {
        public long TotalTransitions
        {
            get; private set;
        }

        /// <summary>
        /// Sets constructor for state machine.
        /// </summary>
        public Func<int, IStateMachine<StateType>> CreateStateMachine
        {
            get; set;
        }

        /// <summary>
        /// Main function. Processes initial states down to final count. 
        /// </summary>
        public virtual BigInteger Process( int bridges, IEnumerable<StateType> initialStates )
        {
            var counts = initialStates.Select( state => new KeyValuePair<StateType, BigInteger>( state, BigInteger.One ) ).ToArray();
            int nn = bridges;
            while (nn > 0)
            {
                counts = Accumulate( CreateStateMachine( --nn ), counts );
            }
            return Total( counts );
        }

        /// <summary>
        /// Detertmines the final total value.
        /// </summary>
        protected BigInteger Total( IEnumerable<KeyValuePair<StateType, BigInteger>> counts )
        {
            BigInteger total = BigInteger.Zero;
            foreach (var kv in counts)
            {
                total += kv.Value;
            }
            return total;
        }

        /// <summary>
        /// Processes one layer.
        /// </summary>
        protected KeyValuePair<StateType, BigInteger>[] Accumulate( IStateMachine<StateType> layer, KeyValuePair<StateType, BigInteger>[] previous )
        {
            var counts = new Dictionary<StateType, BigInteger>();
            long transitions = 0;
            foreach (var kv in previous)
            {
                foreach (var state in layer.Enumerate( kv.Key ))
                {
                    BigInteger n;
                    counts[ state ] = counts.TryGetValue( state, out n ) ? n + kv.Value : kv.Value;
                    ++transitions;
                }
            }
            TotalTransitions += transitions;
            LogProgress( "Layer {0} completed: states={1}, transitions={2}", layer.LayerIndex, counts.Count, transitions );
            return counts.ToArray();
        }

        protected virtual void LogProgress( string format, params object[] args )
        {
        }
    }
    #endregion

    #region DivideAndConquerProcessor
    /// <summary>
    /// Enhanced processor that will partition the work into buckets as required in order to reduce memory load. Works only with basic BigInteger state.
    /// </summary>
    /// <remarks>
    /// Optimal choice of split limit and repartition bucket size is dependent on the number of bridges to be processed as well as memory constraints.
    /// Ideally we want to use as much memory as possible to reduce the number of times the data needs to be repartitioned.
    /// The author has not investigated best practice heuristics.
    /// For better performance we really want to work as much as possible with 64-bit unsigned ints.
    /// (Code for this is really a lot of ugly duplication so is not included here: C# generics does not readily accomodate this.)
    /// This method was successfully used to compute a003516(55) in about 24 hours on a single core, requiring the processing of nearly 400 billion state transitions.
    /// </remarks>
    class DivideAndConquerProcessor : SimpleProcessor<BigInteger>
    {

        public int SplitLimit
        { get; set; }
        public int RepartionBucketSize
        { get; set; }

        /// <summary>
        /// Constructor.
        /// </summary>
        public DivideAndConquerProcessor()
        {
            SplitLimit = 1000000;
            RepartionBucketSize = 150000;
        }

        /// <summary>
        /// Main function. Processes initial states down to final count. 
        /// </summary>
        public override BigInteger Process( int bridges, IEnumerable<BigInteger> initialStates )
        {
            var counts = initialStates.Select( state => new KeyValuePair<BigInteger, BigInteger>( state, BigInteger.One ) ).ToArray();
            return ProcessCluster( new[] { counts }.ToList(), bridges );
        }

        /// <summary>
        /// Called recursively to process a partitioned work-load.
        /// </summary>
        private BigInteger ProcessCluster( List<KeyValuePair<BigInteger, BigInteger>[]> work, int remainingBridges )
        {
            LogProgress( "Begin Cluster at {0}; items={1}", remainingBridges, work.Count );

            BigInteger total = BigInteger.Zero;
            for (int workIndex = 0; workIndex < work.Count; ++workIndex)
            {
                LogProgress( "Layer {0}: Processing List {1} of {2}, items={3}", remainingBridges, workIndex + 1, work.Count, work[ workIndex ].Length );

                var counts = work[ workIndex ];
                work[ workIndex ] = null;

                int b = remainingBridges;
                while (b > 0 && counts.Length < SplitLimit)
                {
                    counts = Accumulate( CreateStateMachine( --b ), counts );
                }
                if (b == 0)
                {
                    total += Total( counts );
                }
                else
                {
                    var partitions = BuildCluster( counts, RepartionBucketSize );
                    GC.Collect();
                    counts = null;
                    total += ProcessCluster( partitions, b );
                }
            }
            LogProgress( "End Cluster at {0}", remainingBridges );
            return total;
        }

        /// <summary>
        /// Given a source list of state-value pairs, partitions the list into multiple buckets for separate processing.
        /// State-value pairs are represented by key-value pairs of BigInteger meander states and NumberType.
        /// </summary>
        private List<KeyValuePair<BigInteger, NumberType>[]> BuildCluster<NumberType>( IEnumerable<KeyValuePair<BigInteger, NumberType>> source, int maxSize )
        {
            var current = new List<List<KeyValuePair<BigInteger, NumberType>>>();
            var output = new List<KeyValuePair<BigInteger, NumberType>[]>();
            current.Add( source.ToList() );
            for (int level = 4; current.Count > 0; ++level)
            {
                var stillToBig = new List<List<KeyValuePair<BigInteger, NumberType>>>();
                for (int currentIndex = 0; currentIndex < current.Count; ++currentIndex)
                {
                    var stateDictionary = current[ currentIndex ];
                    current[ currentIndex ] = null;
                    // divide each item in stateDictionary into a bucket based on MeanderState.Cluster.
                    var index = new Dictionary<int, List<KeyValuePair<BigInteger, NumberType>>>();
                    foreach (var kv in stateDictionary)
                    {
                        var cluster = MeanderState.Cluster( kv.Key, level );
                        List<KeyValuePair<BigInteger, NumberType>> inner;
                        if (!index.TryGetValue( cluster, out inner ))
                        {
                            inner = new List<KeyValuePair<BigInteger, NumberType>>();
                            index.Add( cluster, inner );
                        }
                        inner.Add( kv );
                    }
                    stateDictionary = null;
                    // now sort through each list, putting onto either output list or stillToBig as appropriate.
                    List<KeyValuePair<BigInteger, NumberType>> small = null;
                    foreach (var component in index.Values)
                    {
                        if (component.Count > maxSize)
                        {
                            stillToBig.Add( component );
                        }
                        else if (component.Count > maxSize / 2)
                        {
                            output.Add( component.ToArray() );
                        }
                        else if (small == null)
                        {
                            small = component;
                        }
                        else
                        {
                            small.AddRange( component );
                            if (small.Count > maxSize / 2)
                            {
                                output.Add( small.ToArray() );
                                small = null;
                            }
                        }
                    }
                    if (small != null)
                    {
                        output.Add( small.ToArray() );
                        small = null;
                    }
                }
                current = stillToBig;
            }
            return output;
        }
    }
    #endregion

    #region ArchConfiguration
    /// <summary>
    /// Set of utility methods for exploring arch configurations and meanders.
    /// These are mostly suitable for small n, but have the advantage that the actual meander is obtainable.
    /// </summary>
    /// <remarks>
    /// Both arch configurations and meanders are permutations with a length that is even.
    /// They are represented using integer arrays of their permutation in an obvious manner.
    /// This class although self contained does not include general methods for working with permutations.
    /// </remarks>
    public static class ArchConfiguration
    {

        /// <summary>
        /// The number of arch configurations. (Catalan numbers) (A000108)
        /// </summary>
        public static int ArchConfigurationCount( int n )
        {
            return EnumerateArches( n, arch => 0 ).Count();
        }

        /// <summary>
        /// The number of closed meanders. (A005315)
        /// </summary>
        public static int ClosedMeanderCount( int n )
        {
            return EnumerateArches( n, arch => ArchMeanderCount( arch ) ).Sum();
        }

        /// <summary>
        /// The number of open meanders with an even number of bridges. (A077054)
        /// </summary>
        public static int OpenMeanderCount( int n )
        {
            return EnumerateArches( n, arch => ExteriorArchCount( arch ) * ArchMeanderCount( arch ) ).Sum();
        }

        /// <summary>
        /// The number of symmetric closed meanders. (A060206)
        /// </summary>
        public static int SymmetricClosedMeanderCount( int n )
        {
            int nn = 4 * n - 2;
            return EnumerateArches( nn / 2, arch => arch )
                    .Where( arch => IsSymmetricArch( arch, k => nn - k - 1 ) )
                    .Count();
        }

        /// <summary>
        /// The number of symmetric open meanders. (A223096)
        /// </summary>
        public static int SymmetricOpenMeanderCount( int n )
        {
            int nn = 2 * n + 2;
            return EnumerateArches( nn / 2, arch => arch )
                    .Where( arch => IsSymmetricArch( arch, k => k == 0 ? 0 : nn - k ) )
                    .Count() / 2;
        }

        /// <summary>
        /// The number of arch configurations which have reflective symmetry. (gives A001405)
        /// </summary>
        public static int ReflectiveArchConfigurationCount( int n )
        {
            return EnumerateArches( n, arch => arch )
                .Where( arch => IsReflectiveArch( arch ) )
                .Count();
        }

        /// <summary>
        /// The number of isomorphic arch configurations, where equivalence is defined on rotation. (gives A002995)
        /// </summary>
        public static int CyclicallyDistinctArchConfigurationCount( int n )
        {
            return ArchConfiguratonsByCyclicOrder( n ).Select( kv => kv.Value ).Sum();
        }

        /// <summary>
        /// Groups arch configurations according to cyclic order.
        /// </summary>
        public static SortedDictionary<int, int> ArchConfiguratonsByCyclicOrder( int n )
        {
            var results = EnumerateArches( n, arch => CyclicOrder( arch ) )
                    .GroupBy( k => k ).ToDictionary( g => g.Key, g => g.Count() / g.Key );
            return new SortedDictionary<int, int>( results );
        }

        /// <summary>
        /// Groups arch configurations according to the number of meanders they form. (no obvious use for this break down).
        /// </summary>
        public static SortedDictionary<int, int> ArchConfiguratonsByMeanders( int n )
        {
            var results = EnumerateArches( n, arch => ArchMeanderCount( arch ) )
                    .GroupBy( k => k ).ToDictionary( g => g.Key, g => g.Count() );
            return new SortedDictionary<int, int>( results );
        }

        /// <summary>
        /// Determines the number of meanders an arch can form.
        /// </summary>
        public static int ArchMeanderCount( int[] arch )
        {
            return EnumerateArches( arch.Length / 2, arch2 => arch2 ).Where( arch2 => IsMeander( arch, arch2 ) ).Count();
        }

        /// <summary>
        /// Enumerates all closed meanders with 2n bridges.
        /// </summary>
        public static IEnumerable<T> EnumerateClosedMeanders<T>( int n, Func<int[], int[], T> capture )
        {
            var q = from arch1 in EnumerateArches( n, arch => arch )
                    from arch2 in EnumerateArches( n, arch => arch )
                    where IsMeander( arch1, arch2 )
                    select capture( arch1, arch2 );
            return q;
        }

        /// <summary>
        /// Tests if two arch configurations come together to make a meander.
        /// </summary>
        public static bool IsMeander( int[] arch1, int[] arch2 )
        {
            int n = arch1.Length;
            int p = 0;
            for (int i = 0; i < n; i += 2)
            {
                if (p == 0 && i > 0)
                {
                    return false;
                }
                p = arch2[ arch1[ p ] ];
            }
            if (p != 0)
            {
                throw new ApplicationException( "internal error" );
            }
            return true;
        }

        /// <summary>
        /// Converts lower and upper arch configurations into a meander permutation. Also works with meander systems.
        /// </summary>
        public static int[] ArchesToMeander( int[] arch1, int[] arch2 )
        {
            int n = arch1.Length;
            int[] meander = new int[ n ];
            for (int i = 0; i < n; i += 2)
            {
                meander[ i ] = arch1[ i ];
                meander[ i + 1 ] = arch2[ i + 1 ];
            }
            return meander;
        }

        /// <summary>
        /// Converts a meander into a pair of arches. This is the converse of ArchesToMeander.
        /// </summary>
        public static Tuple<int[], int[]> MeanderToArches( int[] meander )
        {
            int n = meander.Length;
            int[] arch1 = new int[ n ];
            int[] arch2 = new int[ n ];
            for (int i = 0; i < n; i += 2)
            {
                int q = meander[ i ];
                arch1[ i ] = q;
                arch1[ q ] = i;
                int j = i + 1;
                q = meander[ j ];
                arch2[ j ] = q;
                arch2[ q ] = j;
            }

            return Tuple.Create( arch1, arch2 );
        }

        /// <summary>
        /// Returns the number of exterior arches in an arch configuration.
        /// </summary>
        public static int ExteriorArchCount( int[] arch )
        {
            int count = 0;
            for (int p = 0; p < arch.Length; p = arch[ p ] + 1)
            {
                ++count;
            }
            return count;
        }

        /// <summary>
        /// Tests if an arch configuration can be combined with it-self to make a meander.
        /// </summary>
        public static bool IsSymmetricArch( int[] arch, Func<int, int> reflection )
        {
            int n = arch.Length;
            int p = 0;
            for (int i = 0; i < n; i += 2)
            {
                if (p == 0 && i > 0)
                {
                    return false;
                }
                p = arch[ p ];
                p = reflection( arch[ reflection( p ) ] );
            }
            if (p != 0)
            {
                throw new ApplicationException( "internal error" );
            }
            return true;
        }

        /// <summary>
        /// Tests if an arch configuration when reflected gives it-self.
        /// </summary>
        public static bool IsReflectiveArch( int[] arch )
        {
            int nm1 = arch.Length - 1;
            return Enumerable.Range( 0, arch.Length ).All( i => nm1 - arch[ i ] == arch[ nm1 - i ] );
        }

        /// <summary>
        /// Determines the minimum amount of rotation required to bring an arch configuration back to it-self.
        /// </summary>
        public static int CyclicOrder( int[] arch )
        {
            for (int i = 1; i < arch.Length; ++i)
            {
                if (CyclicSelfCompare( arch, i ) == 0)
                {
                    return i;
                }
            }
            return arch.Length;
        }

        /// <summary>
        /// Compares an arch configuration with a cyclic shift of itself.
        /// </summary>
        public static int CyclicSelfCompare( int[] arch, int shift )
        {
            Func<int, int> Cycle = k => (k >= shift ? k : arch.Length + k) - shift;
            int d = 0;
            for (int i = 0; d == 0 && i < arch.Length; ++i)
            {
                d = arch[ Cycle( i ) ].CompareTo( Cycle( arch[ i ] ) );
            }
            return d;
        }

        /// <summary>
        /// Tests if an array of values represents a valid arch configuration permutation.
        /// </summary>
        public static bool IsArchConfiguration( int[] arch )
        {
            for (int i = 0; i < arch.Length; ++i)
            {
                int t = arch[ i ];
                if (t < 0 || t >= arch.Length || t == i || arch[ t ] != i)
                {
                    return false;
                }
            }
            return true;
        }

        /// <summary>
        /// Enumerates all arch configurations with 2n points.
        /// </summary>
        /// <remarks>
        /// This method non-recursively enumerates permuatations which are arch configurations. (Count gives Catalan numbers)
        /// Use this in preference to low level InitialArchConfiguration / NextArchConfiguration methods.
        /// </remarks>
        public static IEnumerable<T> EnumerateArches<T>( int n, Func<int[], T> capture )
        {
            var perm = InitialArchConfiguration( n );
            do
            {
                yield return capture( perm );
            } while (NextArchConfiguration( perm ));
        }

        /// <summary>
        /// Gets an initial arch configuration permutation that can be used with NextArchConfiguration to get the full sequence.
        /// </summary>
        public static int[] InitialArchConfiguration( int n )
        {
            var perm = new int[ 2 * n ];
            for (int i = 0; i < perm.Length; ++i)
            {
                perm[ i ] = perm.Length - 1 - i;
            }
            return perm;
        }

        /// <summary>
        /// Given an arch configuration, updates it to the 'next' one. No state is required to perform this update.
        /// Note that this performs a destructive update for efficiency.
        /// Arch configurations are ordered with the deepest possible nesting first. '((()))' -> '(()())' ... -> '()()()'.
        /// </summary>
        public static bool NextArchConfiguration( int[] perm )
        {
            int pushed = 0;
            int stack = -1;
            int idx = perm.Length;
            bool backtracking = true;
            while (backtracking && idx != 0)
            {
                --idx;
                if (idx == stack)
                {
                    stack = perm[ idx ];
                    --pushed;
                    backtracking = pushed == 0;
                }
                else
                {
                    perm[ perm[ idx ] ] = stack;
                    stack = perm[ idx ];
                    ++pushed;
                }
            }
            if (pushed > 0)
            {
                perm[ idx ] = stack;
                stack = perm[ stack ];
                perm[ perm[ idx ] ] = idx;
                --pushed;
                ++idx;
            }

            do
            {
                if (idx < perm.Length - pushed)
                {
                    perm[ idx ] = stack;
                    stack = idx;
                    ++pushed;
                }
                else
                {
                    perm[ idx ] = stack;
                    stack = perm[ stack ];
                    perm[ perm[ idx ] ] = idx;
                    --pushed;
                }
            } while (++idx < perm.Length);

            return !backtracking;
        }

    }
    #endregion
}