Binomial Coefficient part 2, Catalan Numbers

// Catalan(1.000.000) in 810 ms (Athlon X4 640, XP, 2GB). 
// Amsterdam University Library thanks, there I read
// "Computing Binomial Coefficients" (April 1987!) by P. Goetgheluck.
// The blog about binomials (december 2010) has to be updated.
// "BNM(uint n, uint k)" chooses the fastest algorithm, BNM1 or BNM3. 
// If k is small compared to n, BNM1 will be used.
// Example: BNM1(12,4)=12!/8!/4!=9*10*11*12/(2*3*4)=9*5*11*6/(3*2)
// Otherwise BNM3 will be used, it uses prime power factorization.
// A sieve of Eratosthenes "getComposites" returns a BitArray,
// and an approximation for the square root of n.
// The exact square root is found by "root".
// Goetgheluck's algorithm is used in "getBnmPrimes",
// to get the power "exp(n,k,p)" of a prime.
// A signed integer version "exp(int n, int k, int p)" isn't used,
// but ..... is there for it's own sake. 
// Some care has been taken to work with arrays of uints, ulongs etc as long as possible.
// Result: BNM(6.400.000 , 2.133.333) in 3510 ms, that's nearly 10 % faster.

using Xint = System.Numerics.BigInteger;
using System.Collections.Generic;
using System.Threading.Tasks;
using System.Collections;
using System.Diagnostics;
using System;
class Binomial
{

    public static Xint CAT(uint n)
    {
        return n < 3 ? n / 2 + 1 : n < 11 ? BNM1(n * 2, n++) / n : BNM3(n * 2, n++) / n;
    }

    public static Xint BNM(uint n, uint k)
    {
        if (k > n) return 0;
        if (k > n / 2) k = n - k;
        if (k == 0) return 1;
        if (k == 1) return n;
        if (k == 2) return n <= ushort.MaxValue ? n-- * n / 2 : (ulong)(n--) * n / 2;
        if (k < 11) return BNM1(n, k);
        if (n <= 64) return BNM3(n, k);
        if (n <= 128) return k < 12 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 256) return k < 12 + 5 * (n - 128) / 128 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 512) return k < 17 + 2 * (n - 256) / 256 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 1024) return k < 19 + 9 * (n - 512) / 512 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 2048) return k < 28 + 13 * (n - 1024) / 1024 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 4096) return k < 41 + 17 * (n - 2048) / 2048 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 8192) return k < 58 + 31 * (n - 4096) / 4096 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 16384) return k < 89 + 45 * (n - 8192) / 8192 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 32768) return k < 134 + 60 * (n - 16384) / 16384 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 65536) return k < 194 + 88 * (n - 32768) / 32768 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 131072) return k < 282 + 120 * (n - 65536) / 65536 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 262144) return k < 402 + 162 * (n - 131072) / 131072 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 524288) return k < 564 + 229 * (n - 262144) / 262144 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 1048576) return k < 793 + 312 * (n - 524288) / 524288 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 2097152) return k < 1105 + 430 * (n - 1048576) / 1048576 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 4194304) return k < 1535 + 583 * (n - 2097152) / 2097152 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 8388608) return k < 2118 + 806 * (n - 4194304) / 4194304 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 16777216) return k < 2924 + (ulong)(1700) * (n - 8388608) / 8388608 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 33554432) return k < 4624 + (ulong)(2328) * (n - 16777216) / 16777216 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 67108864) return k < 6952 + (ulong)(3135) * (n - 33554432) / 33554432 ? BNM1(n, k) : BNM3(n, k);
        if (n <= 134217728) return k < 10087 + (ulong)(4282) * (n - 67108864) / 67108864 ? BNM1(n, k) : BNM3(n, k);
        return k < 10087 + 4282 * (n >> 13) / 16384 ? BNM1(n, k) : BNM3(n, k);
    }

    private static Xint BNM1(uint n, uint k)
    {
        ulong uu = 0;
        Xint U = 0;
        uint i = n - k + 1;
        uint u = i;
        if ((i++ & 1) == 0)
        {
            u /= 2;
            goto L1;
        }
    L0: if (i <= n && u <= ushort.MaxValue && i <= ushort.MaxValue)
        {
            u *= (i++ / 2);
        }
        else
        {
            if (i > n) goto M0;
            else
            {
                uu = (ulong)(u) * (i++ / 2);
                goto L2;
            }
        }
    L1: if (i <= n && u <= ushort.MaxValue && i <= ushort.MaxValue)
        {
            u *= i++;
            goto L0;
        }
        else
        {
            if (i > n) goto M0;
            else
            {
                uu = (ulong)(u) * i++;
                goto L3;
            }
        }

    L2: if (i <= n && uu <= uint.MaxValue)
        {
            uu *= i++;
        }
        else
        {
            if (i > n) goto M0;
            else
            {
                U = (Xint)(uu) * i++;
                goto L4;
            }
        }
    L3: if (i <= n && uu <= uint.MaxValue)
        {
            uu *= (i++ / 2);
            goto L2;
        }
        else
        {
            if (i > n) goto M0;
            else
            {
                U = (Xint)(uu) * (i++ / 2);
                goto L5;
            }
        }

    L4: if (i <= n)
        {
            U *= (i++ / 2);
        }
        else goto M0;
    L5: if (i <= n)
        {
            U *= i++;
            goto L4;
        }

    M0: switch (k)
        {
            case 3:
                if (uu == 0) return u / 3 << 1 - (int)(n & 1);
                if (U.IsZero) return uu / 3 << 1 - (int)(n & 1);
                return U / 3 << 1 - (int)(n & 1);
            case 4:
                if (uu == 0) return u / 6;
                if (U.IsZero) return uu / 6;
                return U / 6;
            case 5:
                if (uu == 0) return (u >> (int)(n & 1)) / 15;
                if (U.IsZero) return (uu >> (int)(n & 1)) / 15;
                return (U >> (int)(n & 1)) / 15;
            case 6:
                if (uu == 0) return u / 90;
                if (U.IsZero) return uu / 90;
                return U / 90;
            case 7:
                if (U.IsZero) return (uu >> (int)(n & 1)) / 315;
                return (U >> (int)(n & 1)) / 315;
            case 8:
                if (U.IsZero) return uu / 2520;
                return U / 2520;
            case 9:
                if (U.IsZero) return (uu >> (int)(n & 1)) / 11340;
                return (U >> (int)(n & 1)) / 11340;
            case 10:
                if (U.IsZero) return uu / 113400;
                return U / 113400;
            case 11:
                return (U >> (int)(n & 1)) / 623700;
            case 12:
                return U / 7484400;
            case 13:
                return (U >> (int)(n & 1)) / 48648600;
            case 14:
                return U / 681080400;
            case 15:
                return (U >> 3 + (int)(n & 1)) / 638512875;
            case 16:
                return (U >> 7) / 638512875;
            case 17:
                return (U >> 6 + (int)(n & 1)) / 10854718875;
            case 18:
                return (U >> 7) / 97692469875;
            case 19:
                return (U >> 6 + (int)(n & 1)) / 1856156927625;
            case 20:
                return (U >> 8) / 9280784638125;
            case 21:
                return (U >> 7 + (int)(n & 1)) / 194896477400625;
            case 22:
                return (U >> 8) / 2143861251406875;
            case 23:
                return (U >> 7 + (int)(n & 1)) / 49308808782358125;
            case 24:
                return (U >> 10) / 147926426347074375;
            case 25:
                return (U >> 9 + (int)(n & 1)) / 3698160658676859375;
            default:
                Xint V = 3698160658676859375;
                i = 26;
            M1: if (i <= k) V *= (i++ / 2);
                else return (U >> 9 + (int)(n & 1)) / V;
                if (i <= k)
                {
                    V *= i++;
                    goto M1;
                }
                return (U >> 10) / V;
        }
    }

    private static Xint BNM3(uint n, uint k)
    {
        if (n <= 65536)
        {
            uint[] f = getBnmPrimes(n, k);
            int i = f.Length;
            if (i == 1) return (Xint)(f[0]) << exp(n, k, 2);
            int j = 1 << fL2(i);
            uint max = 0;
            if (i != j)
            {
                f = uintSpecialProduct(f, i, j, out max);
            }
            for (; j > 1 && max <= ushort.MaxValue; j /= 2) f = uintProduct(f, j, out max);
            if (j == 1) return (Xint)(f[0]) << exp(n, k, 2);
            ulong mmax = 0;
            ulong[] ff = ulongProduct(f, j, out mmax);
            j /= 2;
            if (j == 1) return (Xint)(ff[0]) << exp(n, k, 2);
            Array.Resize(ref f, 0);
            for (; j > 1 && mmax <= uint.MaxValue; j /= 2) ff = ulongProduct(ff, j, out mmax);
            if (j == 1) return (Xint)(ff[0]) << exp(n, k, 2);
            Xint[] F = XintProduct(ff, j);
            j /= 2;
            if (j == 1) return F[0] << exp(n, k, 2);
            Array.Resize(ref ff, 0);
            for (; j > 1 && F[0].ToByteArray().Length < 376; j /= 2) F = SmallProduct(F, j);
            for (; j > 1; j /= 2) F = LargeProduct(F, j);
            return F[0] << exp(n, k, 2);
        }
        else
        {
            uint[] f = getBnmPrimes(n, k);
            int i = f.Length;
            if (i == 1) return (Xint)(f[0]) << exp(n, k, 2);
            int j = 1 << fL2(i);
            ulong mmax = 0;
            ulong[] ff;
            if (i != j)
            {
                ff = ulongSpecialProduct(f, i, j, out mmax);
            }
            else
            {
                ff = ulongProduct(f, j, out mmax);
                j /= 2;
            }
            if (j == 1) return (Xint)(ff[0]) << exp(n, k, 2);
            Array.Resize(ref f, 0);
            for (; j > 1 && mmax <= uint.MaxValue; j /= 2) ff = ulongProduct(ff, j, out mmax);
            if (j == 1) return (Xint)(ff[0]) << exp(n, k, 2);
            Xint[] F = XintProduct(ff, j);
            j /= 2;
            if (j == 1) return F[0] << exp(n, k, 2);
            Array.Resize(ref ff, 0);
            for (; j > 1 && F[0].ToByteArray().Length < 376; j /= 2) F = SmallProduct(F, j);
            for (; j > 1; j /= 2) F = LargeProduct(F, j);
            return F[0] << exp(n, k, 2);
        }
    }

    private static uint[] getBnmPrimes(uint n, uint k)
    {
        int i;
        BitArray composite = getComposites(n, out i);
        uint rt = root(n, i);
        List<uint> primes = new List<uint>();
        int e = exp(n, k, 3);
        for (i = 0; i < e; i++) primes.Add(3);

        uint p = 5;
        i = 0;
        int j = 0;

    L0: if (p > rt) { e = (int)(n / 2); goto L1; }
        if (!composite[i++])
        {
            e = exp(n, k, p);
            for (j = 0; j < e; j++) primes.Add(p);
        }
        p += 2;
        if (p > rt) { e = (int)(n / 2); goto L2; }
        if (!composite[i++])
        {
            e = exp(n, k, p);
            for (j = 0; j < e; j++) primes.Add(p);
        }
        p += 4;
        goto L0;

    L1: if (p > e) goto L3;
        if (!composite[i++] && n % p < k % p) primes.Add(p);
        p += 2;
    L2: if (p > e) goto L3;
        if (!composite[i++] && n % p < k % p) primes.Add(p);
        p += 4;
        goto L1;

    L3: p = n - k + 1;
        p += 1 - (p & 1);
        p += p % 6 & 2;           // if (p % 6 == 3) p += 2;
        i = (int)(p / 3 - 1);
        if ((i & 1) != 0) goto L5;

    L4: if (p > n) goto L6;
        if (!composite[i++]) primes.Add(p);
        p += 2;
    L5: if (p > n) goto L6;
        if (!composite[i++]) primes.Add(p);
        p += 4;
        goto L4;

    L6: return primes.ToArray();
    }

    private static BitArray getComposites(uint n, out int i)
    {
        int len = (int)(n / 3);
        BitArray composite = new BitArray(len);
        i = 0;
        int d1 = 8, d2 = 0, p1 = 3, p2 = 1, s = 7, inc = 10;
        while (s < len)
        {
            if (!composite[i++])
            {
                int j = s, m = s + p1;
                for (; m < len; j += inc, m += inc)
                {
                    composite[j] = true;
                    composite[m] = true;
                }
                if (j < len)
                    composite[j] = true;
            }
            d2 += 8; s += d2; p2 += 8; inc += 4;
            if (!composite[i++])
            {
                int j = s, m = s + p2;
                for (; m < len; j += inc, m += inc)
                {
                    composite[j] = true;
                    composite[m] = true;
                }
                if (j < len)
                    composite[j] = true;
            }
            d1 += 16; s += d1; p1 += 4; inc += 8;
        }
        return composite;
    }

    private static uint root(uint n, int i)
    {
        uint rt = (uint)(i * 3 + 3);
        uint sq = rt * rt;
        for (; sq < n; sq += 2 * rt | 1, rt++) ;
        for (; sq > n; rt--, sq -= 2 * rt | 1) ;
        return rt;
    }

    private static int exp(uint n, uint k, uint p)
    {
        int e = 0;

    L0: uint qn = n / p;
        uint qk = k / p;
        if ((qn - qk) * p > n - k)
        {
            e++;
            n = qn;
            if (qk == 0) goto L2;
            k = qk;
            goto L1;
        }
        if (qk == 0) return e;
        n = qn;
        k = qk;
        goto L0;

    L1: qn = n / p;
        qk = k / p;
        if ((qn - qk) * p >= n - k)
        {
            e++;
            n = qn;
            if (qk == 0) goto L2;
            k = qk;
            goto L1;
        }
        if (qk == 0) return e;
        n = qn;
        k = qk;
        goto L0;

    L2: qn = n / p;
        if (n == qn * p) e++;
        else return e;
        n = qn;
        goto L2;
    }

    private static int exp(int n, int k, int p)
    {
        int e = 0, rn = 0, rk = 0;

    L0: n = Math.DivRem(n, p, out rn);
        k = Math.DivRem(k, p, out rk);
        if (rk > rn)
        {
            e++;
            if (k == 0) goto L2;
            goto L1;
        }
        if (k == 0) return e;
        goto L0;

    L1: n = Math.DivRem(n, p, out rn);
        k = Math.DivRem(k, p, out rk);
        if (rk >= rn)
        {
            e++;
            if (k == 0) goto L2;
            goto L1;
        }
        if (k == 0) return e;
        goto L0;

    L2: n = Math.DivRem(n, p, out rn);
        if (rn == 0) e++;
        else return e;
        goto L2;
    }

    private static uint[] uintSpecialProduct(uint[] f, int i, int j, out uint max)
    {
        max = 0;
        uint p = 0;
        uint[] newF = new uint[j];
        int m = i - j;                                // ?!?!
        i = j - m;                                    // ?!?!
        int n = 0;
        while (n < i) newF[n++] = f[m++];
        i = 0;
        while (n < j)
        {
            p = f[m++] * f[i++];
            if (max < p) max = p;
            newF[n++] = p;
        }
        return newF;
    }
    private static ulong[] ulongSpecialProduct(uint[] f, int i, int j, out ulong mmax)
    {
        mmax = 0;
        ulong pp = 0;
        ulong[] newF = new ulong[j];
        int m = i - j;
        i = j - m;
        int n = 0;
        while (n < i) newF[n++] = f[m++];
        i = 0;
        while (n < j)
        {
            pp = (ulong)(f[m++]) * f[i++];
            if (mmax < pp) mmax = pp;
            newF[n++] = pp;
        }
        return newF;
    }
    private static uint[] uintProduct(uint[] f, int j, out uint max)
    {
        max = 0;
        uint p = 0;
        int k = j-- / 2;
        uint[] newF = new uint[k];
        for (int i = 0; i < k; i++, j--)
        {
            p = f[i] * f[j];
            if (max < p) max = p;
            newF[i] = p;
        }
        return newF;
    }
    private static ulong[] ulongProduct(uint[] f, int j, out ulong mmax)
    {
        mmax = 0;
        ulong pp = 0;
        int k = j-- / 2;
        ulong[] newF = new ulong[k];
        for (int i = 0; i < k; i++, j--)
        {
            pp = (ulong)(f[i]) * f[j];
            if (mmax < pp) mmax = pp;
            newF[i] = pp;
        }
        return newF;
    }
    private static ulong[] ulongProduct(ulong[] ff, int j, out ulong mmax)
    {
        mmax = 0;
        ulong pp = 0;
        int k = j-- / 2;
        ulong[] newF = new ulong[k];
        for (int i = 0; i < k; i++, j--)
        {
            pp = ff[i] * ff[j];
            if (mmax < pp) mmax = pp;
            newF[i] = pp;
        }
        return newF;
    }
    private static Xint[] XintProduct(ulong[] ff, int j)
    {
        int k = j-- / 2;
        Xint[] NewF = new Xint[k];
        for (int i = 0; i < k; i++, j--) NewF[i] = (Xint)(ff[i]) * ff[j];
        return NewF;
    }
    private static Xint[] SmallProduct(Xint[] F, int j)
    {
        int k = j-- / 2;
        Xint[] NewF = new Xint[k];
        for (int i = 0; i < k; i++, j--) NewF[i] = F[i] * F[j];
        return NewF;
    }
    private static Xint[] LargeProduct(Xint[] F, int j)
    {
        int k = j-- / 2;
        Xint[] NewF = new Xint[k];
        for (int i = 0; i < k; i++, j--) NewF[i] = MTP(F[i], F[j]);
        return NewF;
    }

    public static int bL(Xint U)
    {
        byte[] bytes = (U.Sign * U).ToByteArray();
        int i = bytes.Length - 1;
        return i * 8 | bitLengthMostSignificantByte(bytes[i]);
    }
    private static int bitLengthMostSignificantByte(byte b)
    {
        return b < 08 ? b < 02 ? b < 01 ? 0 : 1 :
                                 b < 04 ? 2 : 3 :
                        b < 32 ? b < 16 ? 4 : 5 :
                                 b < 64 ? 6 : 7;
    }

    private static Xint MTP(Xint U, Xint V)
    {
        return MTP(U, V, Xint.Max(U.Sign * U, V.Sign * V).ToByteArray().Length << 3);
    }
    private static Xint MTP(Xint U, Xint V, int n)
    {
        if (n <= 3000) return U * V;
        if (n <= 6000) return TC2(U, V, n);
        if (n <= 10000) return TC3(U, V, n);
        if (n <= 40000) return TC4(U, V, n);
        return TC2P(U, V, n);
    }
    private static Xint MTPr(Xint U, Xint V, int n)
    {
        if (n <= 3000) return U * V;
        if (n <= 6000) return TC2(U, V, n);
        if (n <= 10000) return TC3(U, V, n);
        return TC4(U, V, n);
    }
    private static Xint TC2(Xint U1, Xint V1, int n)
    {
        n >>= 1;
        Xint Mask = (Xint.One << n) - 1;
        Xint U0 = U1 & Mask; U1 >>= n;
        Xint V0 = V1 & Mask; V1 >>= n;
        Xint P0 = MTPr(U0, V0, n);
        Xint P2 = MTPr(U1, V1, n);
        return ((P2 << n) + (MTPr(U0 + U1, V0 + V1, n) - (P0 + P2)) << n) + P0;
    }
    private static Xint TC3(Xint U2, Xint V2, int n)
    {
        n = (int)((long)(n) * 0x55555556 >> 32); // n /= 3;
        Xint Mask = (Xint.One << n) - 1;
        Xint U0 = U2 & Mask; U2 >>= n;
        Xint U1 = U2 & Mask; U2 >>= n;
        Xint V0 = V2 & Mask; V2 >>= n;
        Xint V1 = V2 & Mask; V2 >>= n;
        Xint W0 = MTPr(U0, V0, n);
        Xint W4 = MTPr(U2, V2, n);
        Xint P3 = MTPr((((U2 << 1) + U1) << 1) + U0, (((V2 << 1) + V1 << 1)) + V0, n);
        U2 += U0;
        V2 += V0;
        Xint P2 = MTPr(U2 - U1, V2 - V1, n);
        Xint P1 = MTPr(U2 + U1, V2 + V1, n);
        Xint W2 = (P1 + P2 >> 1) - (W0 + W4);
        Xint W3 = W0 - P1;
        W3 = ((W3 + P3 - P2 >> 1) + W3) / 3 - (W4 << 1);
        Xint W1 = P1 - (W4 + W3 + W2 + W0);
        return ((((W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0;
    }
    private static Xint TC4(Xint U3, Xint V3, int n)
    {
        n >>= 2;
        Xint Mask = (Xint.One << n) - 1;
        Xint U0 = U3 & Mask; U3 >>= n;
        Xint U1 = U3 & Mask; U3 >>= n;
        Xint U2 = U3 & Mask; U3 >>= n;
        Xint V0 = V3 & Mask; V3 >>= n;
        Xint V1 = V3 & Mask; V3 >>= n;
        Xint V2 = V3 & Mask; V3 >>= n;

        Xint W0 = MTPr(U0, V0, n);                               //  0
        U0 += U2; U1 += U3;
        V0 += V2; V1 += V3;
        Xint P1 = MTPr(U0 + U1, V0 + V1, n);                     //  1
        Xint P2 = MTPr(U0 - U1, V0 - V1, n);                     // -1
        U0 += 3 * U2; U1 += 3 * U3;
        V0 += 3 * V2; V1 += 3 * V3;
        Xint P3 = MTPr(U0 + (U1 << 1), V0 + (V1 << 1), n);       //  2
        Xint P4 = MTPr(U0 - (U1 << 1), V0 - (V1 << 1), n);       // -2
        Xint P5 = MTPr(U0 + 12 * U2 + ((U1 + 12 * U3) << 2),
                       V0 + 12 * V2 + ((V1 + 12 * V3) << 2), n); //  4
        Xint W6 = MTPr(U3, V3, n);                               //  inf

        Xint W1 = P1 + P2;
        Xint W4 = (((((P3 + P4) >> 1) - (W1 << 1)) / 3 + W0) >> 2) - 5 * W6;
        Xint W2 = (W1 >> 1) - (W6 + W4 + W0);
        P1 = P1 - P2;
        P4 = P4 - P3;
        Xint W5 = ((P1 >> 1) + (5 * P4 + P5 - W0 >> 4) - ((((W6 << 4) + W4) << 4) + W2)) / 45;
        W1 = ((P4 >> 2) + (P1 << 1)) / 3 + (W5 << 2);
        Xint W3 = (P1 >> 1) - (W1 + W5);
        return ((((((W6 << n) + W5 << n) + W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0;
    }
    private static Xint TC2P(Xint A, Xint B, int n)
    {
        n >>= 1;
        Xint Mask = (Xint.One << n) - 1;
        Xint[] U = new Xint[3];
        U[0] = A & Mask; A >>= n; U[2] = A; U[1] = U[0] + A;
        Xint[] V = new Xint[3];
        V[0] = B & Mask; B >>= n; V[2] = B; V[1] = V[0] + B;
        Xint[] P = new Xint[3];
        Parallel.For(0, 3, (int i) => P[i] = MTPr(U[i], V[i], n));
        return ((P[2] << n) + P[1] - (P[0] + P[2]) << n) + P[0];
    }

    private static int fL2(int i)
    {
        return
        i < 1 << 15 ? i < 1 << 07 ? i < 1 << 03 ? i < 1 << 01 ? i < 1 << 00 ? -1 : 00 :
                                                                i < 1 << 02 ? 01 : 02 :
                                                  i < 1 << 05 ? i < 1 << 04 ? 03 : 04 :
                                                                i < 1 << 06 ? 05 : 06 :
                                    i < 1 << 11 ? i < 1 << 09 ? i < 1 << 08 ? 07 : 08 :
                                                                i < 1 << 10 ? 09 : 10 :
                                                  i < 1 << 13 ? i < 1 << 12 ? 11 : 12 :
                                                                i < 1 << 14 ? 13 : 14 :
                      i < 1 << 23 ? i < 1 << 19 ? i < 1 << 17 ? i < 1 << 16 ? 15 : 16 :
                                                                i < 1 << 18 ? 17 : 18 :
                                                  i < 1 << 21 ? i < 1 << 20 ? 19 : 20 :
                                                                i < 1 << 22 ? 21 : 22 :
                                    i < 1 << 27 ? i < 1 << 25 ? i < 1 << 24 ? 23 : 24 :
                                                                i < 1 << 26 ? 25 : 26 :
                                                  i < 1 << 29 ? i < 1 << 28 ? 27 : 28 :
                                                                i < 1 << 30 ? 29 : 30;
    }

    private static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        BNM(1000000, 500000);
        sw.Restart();
        BNM(6400000, 2133333);
        sw.Stop();
        Console.WriteLine("BNM(6.400.000 , 2.133.333)");
        Console.WriteLine(sw.ElapsedMilliseconds + " ms");
        Console.WriteLine();

        sw.Restart();
        Xint C = CAT(100);
        sw.Stop();
        Console.WriteLine("Catalan(100)");
        Console.WriteLine(sw.ElapsedMilliseconds + " ms");
        Console.WriteLine(bL(C) + " bits");
        Console.WriteLine(C);
        Console.WriteLine();

        sw.Restart();
        C = CAT(100);
        sw.Stop();
        Console.WriteLine("Catalan(100)");
        Console.WriteLine(sw.ElapsedMilliseconds + " ms");
        Console.WriteLine(bL(C) + " bits");
        Console.WriteLine(C);
        Console.WriteLine();

        sw.Restart();
        C = CAT(200);
        sw.Stop();
        Console.WriteLine("Catalan(200)");
        Console.WriteLine(sw.ElapsedMilliseconds + " ms");
        Console.WriteLine(bL(C) + " bits");
        Console.WriteLine(C);
        Console.WriteLine();

        sw.Restart();
        for (int i = 0; i < 1000; i++)
        {
            C = CAT(1000);
        }
        sw.Stop();
        Console.WriteLine("Catalan(1000)");
        Console.WriteLine(sw.ElapsedMilliseconds + " us");
        Console.WriteLine(bL(C) + " bits");
        Console.WriteLine(C);
        Console.WriteLine();

        sw.Restart();
        C = CAT(1000000);
        sw.Stop();
        Console.WriteLine("Catalan(1.000.000)");
        Console.WriteLine(sw.ElapsedMilliseconds + " ms");
        Console.WriteLine(bL(C) + " bits");
        Console.WriteLine();

        sw.Restart();
        C = CAT(2000000);
        sw.Stop();
        Console.WriteLine("Catalan(2.000.000)");
        Console.WriteLine(sw.ElapsedMilliseconds + " ms");
        Console.WriteLine(bL(C) + " bits");
        Console.ReadLine();
    }
}