Pagina's

2012/12/29

A218076 extended with C++

/*
The next term is: F(127495932) with 53 consecutive zeros.
I used Visual C++ 2010 Express, GMP, and MPIR,
see: GMP on windows

The C++ code below is about five times faster than the C# code before.
And... the next term will be larger than 150 * 10^6. 
*/


#include "stdafx.h"
#include "C:\projects\mpir-2.6.0\lib\Win32\Release\gmpxx.h"
#include "ctime"
#include "iostream"
using namespace std;

int lz(unsigned __int8 b)
{
    return b < 1 << 4 ? b < 1 << 2 ? b < 1 << 1 ? 7 : 6 :
                                     b < 1 << 3 ? 5 : 4 :
                        b < 1 << 6 ? b < 1 << 5 ? 3 : 2 :
                                     b < 1 << 7 ? 1 : 0;
}

int tz(unsigned __int8 b)
{
    return 7 - lz((unsigned __int8)(b & -(__int8)b));
}

void main()
{
    int i, i0, i1, j, k, z, zmax, size;
    cout << " low: "cin >> i0;
    cout << "high: "cin >> i1;
    cout << "zmax: "cin >> zmax;
    k = (zmax - 6) / 8;
    mpz_t F0, F1;
    mpz_inits(F0, F1, NULL);
    time_t now = time(0); cout << ctime(&now);
    cout << "compute F(" << dec << i1 << ")" << endl ;
    mpz_fib_ui(F0, i1); // mpz_out_str(stdout,16,F0)
    now = time(0); cout << ctime(&now);
    mpz_set(F1, F0);
    cout << "compute F's(" << dec << i0 - 2 << " & " << dec << i0 - 1 << ")" << endl;
    mpz_fib2_ui(F1, F0, i0 - 1);
    now = time(0); cout << ctime(&now) << endl;
    cout << "enter a number: ";
    cin >> i;
    now = time(0); cout << ctime(&now);
    cout << "S T A R T" << endl << endl;
    unsigned __int8 *p , *q;
    for(i = i0; i <= i1; i++)
    {
        if(i % 100000 == 0
        {
            now = time(0); cout << ctime(&now);
            cout << "F(" << dec << i << ")" << endl << endl;
        }
        mpz_add(F0,F0,F1);
        p = (unsigned __int8 *)F0[0]._mp_d;
        j = 1;
        if (*p == 0)
        {
            for (z = 8, p++, j++; *p == 0; p++, j++, z += 8);
            if (zmax < z + 7)
            {
                z += tz(*p);
                if (zmax < z)
                {
                    zmax = z;
                    k = (z - 6) / 8;
                    now = time(0); cout << "  " << ctime(&now);
                    cout << "  LOWZS: " << dec << z <<"  F(" << dec << i << ")" << endl << endl;
                }
            }
        }
        size = mpz_sizeinbase(F0, 256);
        for (j += k, p += k; j < size; j += k, p += k)
        {
            if (*p == 0)
            {
                for (z = 8, q = p - 1; *q == 0; q--, z += 8);
                for (p++, j++; *p == 0; p++, j++, z += 8);
                if (zmax < z + 14)
                {
                    z += lz(*q);
                    if (zmax < z + 7)
                    {
                        z += tz(*p);
                        if (zmax < z)
                        {
                            zmax = z;
                            k = (z - 6) / 8;
                            now = time(0); cout << "  "  << ctime(&now);
                            cout << "  ZEROS: " << dec << z << "  F(" << dec << i << ") :ZEROS" << endl << endl;
                        }
                    }
                 }
            }
        }
        mpz_swap(F0,F1);
    }
    now = time(0); cout << ctime(&now);
    cout << "F(" << dec << i1 << ")" << endl;
    cout << "enter a number: ";
    cin >> i;
    now = time(0); cout << ctime(&now);
    cout << "compute F(" << dec << i1 << ")" << endl ;
    mpz_fib_ui(F0, i1);
    now = time(0); cout << ctime(&now) << endl;
    if (mpz_cmp(F0, F1) != 0cout << endl << "E R R O R ! ! ! !" << endl;
    mpz_clears(F0, F1, NULL);
    cout << "zmax: " << dec << zmax << endl;
    cout << endl << "R E A D Y";
    cin >> i;
}

2012/11/24

A218076 extended

// The next two terms are:
// F(80141430) 49 consecutive zeros, and
// F(84300971) 51 consecutive zeros
// I used Emil Stefanov's wrapper for GMP 
// (The GNU Multiple Precision Arithmetic Library).
// www.emilstefanov.net/Projects/GnuMpDotNet/
// Copied "libgmp-3.dll" to "C:\WINDOWS\System32\"
// This general compilation computes F(10^9) in 42 seconds.
// Two copies of the program below, for example one
// scanning the range 80000000-82000000, the other 
// scanning the range 82000000-84000000, take ~32 hours.
// And... the next term will be larger than 95999999.

using System;
using System.Threading.Tasks;
using Emil.GMP;  // Project>>Add Reference>>Emil.GMP.dll
class F_cz_gmp
{
    static void Main()
    {
        int i = 80000000;                          // 82000000
        int zmax = 48;
        int[] z = { 0, 0 };
        int k = (zmax - 6) / 8;
        Console.WriteLine("first i: " + i);
        BigInt[] F = new BigInt[2];
        F[1] = BigInt.Fibonacci(i - 1, out F[0]);
        for (; i < 82000000; i += 2)               // 84000000
        {
            if (i % 500000 == 0) Console.WriteLine(i + " " + DateTime.Now);
            F[0] += F[1];
            F[1] += F[0];
            Parallel.For(0, 2, (int j) => z[j] = zeros(F[j], zmax, k));
            GC.Collect();
            if (z[0] > zmax)
            {
                zmax = z[0];
                k = (zmax - 6) / 8;
                Console.WriteLine(i + " " + zmax);
            }
            if (z[1] > zmax)
            {
                zmax = z[1];
                k = (zmax - 6) / 8;
                Console.WriteLine((i + 1) + " " + zmax);
            }
        }
        Console.WriteLine("last i: " + (i - 1) + " " + DateTime.Now);
        Console.ReadLine();
    }
    private static int zeros(BigInt F, int zmax, int k)
    {
        byte[] b = F.ToByteArray(); // F.ToUintArray ~4 times faster
        int bL = b.Length;
        int j = 0;
        int z;
        if (b[0] == 0)
        {
            z = 8;
            while (b[++j] == 0) z += 8; // add zerobits
            if (zmax < z + 7)
            {
                z += tz(b[j]);
                if (zmax < z)
                {
                    zmax = z;
                    k = (z - 6) / 8;
                }
            }
        }
        j += k;
        while (j < bL)
        {
            if (b[j] == 0)
            {
                z = 8;
                int j0 = j;
                while (b[--j0] == 0) z += 8; // add zerobits
                while (b[++j] == 0) z += 8;
                if (zmax < z + 14)
                {
                    z += lz(b[j0]);
                    if (zmax < z + 7)
                    {
                        z += tz(b[j]);
                        if (zmax < z)
                        {
                            zmax = z;
                            k = (z - 6) / 8;
                        }
                    }
                }
            }
            j += k;
        }
        return zmax;
    }
    private static int lz(byte b) // leading zeros, b > 0
    {
        return b < 1 << 4 ? b < 1 << 2 ? b < 1 << 1 ? 7 : 6 :
                                         b < 1 << 3 ? 5 : 4 :
                            b < 1 << 6 ? b < 1 << 5 ? 3 : 2 :
                                         b < 1 << 7 ? 1 : 0;
    }
    private static int tz(byte b) // trailing zeros
    {
        return 7 - lz((byte)(b & -(sbyte)b));
    }
}

2012/10/31

Consecutive zeros in Fibonacci numbers

// When I looked at the Fibonacci Sequence Binary Plot
// http://mathworld.wolfram.com/FibonacciNumber.html
// and the statement: F(2^n+2^(n+1)) ends in n+2 zeros.
// I wondered: Are there n+2 consecutive zeros for smaller indexes?
// Or: Indexes of Fibonacci numbers whose binary expansions have 
//     record numbers of consecutive zeros.
// It seems there are.
// Example: The first four records occur at 3, 6, 12, and 19:
// F(3)=10 2 (one zero)
// F(6)=1000 2 (three zeros)
// F(12)=10010000 2 (four zeros)
// F(19)=1000001010101 2 (five zeros)
//
//      Index  Zeros 
//          3    1  
//          6    3  
//         12    4  
//         19    5
//         38    6
//         42    7
//         68   12
//        243   13
//        384   14  
//        515   16
//        740   19
//       1709   24
//       5151   27
//      11049   28
//      45641   30
//      94729   31
//     185610   34
//     644593   35
//     726681   39
//    2296396   40
//    3098358   42
//    6178778   43
//   15743325   45
//   22436908   48
//
// In the program below, 
// the first loop "for (; zmax < 15; i++)" 
// searches for consecutive zerobits, bit by bit.
// As soon as 15 zeros are found, 
// the second loop "for (; i < int.MaxValue; i++)"
// searches for consecutive zerobytes, byte by byte,
// adding leading/trailing zerobits of the bytes before/after those zerobytes.
// It's much faster than searching bit by bit.

using System;
using Xint = System.Numerics.BigInteger; // Project>>Add Reference>>System.Numerics
class fibonacciRepeatingZeros
{
    static void Main()
    {
        Xint F0 = 1, F1 = 0, F2 = 0;
        int i = 1, z, zi = 0, zmax = 0; // z=zerobits,zi=maxzerobitsF(i)
        for (; zmax < 15; i++)
        {
            F2 = F1 + F0; F0 = F1; F1 = F2;
            z = 0;
            for (; !F2.IsOne; F2 /= 2)
            {
                if (F2.IsEven)
                {
                    if (zi < ++z) zi = z;
                }
                else z = 0;
            }
            if (zmax < zi)
            {
                zmax = zi;
                Console.WriteLine(i + " " + zi);
            }
        }
        byte[] b; int j, j0, bL; z = 0;
        for (; i < int.MaxValue; i++)
        {
            F2 = F1 + F0; F0 = F1; F1 = F2;
            b = F2.ToByteArray();
            bL = b.Length - 1;
            if (b[bL] == 0) bL--;
            j = 0;
            if (b[0] == 0)
            {
                z = 8;
                while (b[++j] == 0) z += 8; // add zerobits
                if (zi < z + 7)
                {
                    z += tz(b[j]);
                    if (zi < z) zi = z;
                }
                z = 0;
            }
            while (++j < bL)
            {
                if (b[j] == 0)
                {
                    j0 = j - 1;
                    z = 8;
                    while (b[++j] == 0) z += 8; // add zerobits
                    if (zi < z + 14)
                    {
                        z += lz(b[j0]);
                        if (zi < z + 7)
                        {
                            z += tz(b[j]);
                            if (zi < z) zi = z;
                        }
                    }
                    z = 0;
                }
            }
            if (zmax < zi)
            {
                zmax = zi;
                Console.WriteLine(i + " " + zi);
            }
        }
        Console.ReadLine();
    }
    private static int lz(byte b) // leading zeros, b > 0
    {
        return b < 1 << 4 ? b < 1 << 2 ? b < 1 << 1 ? 7 : 6 :
                                         b < 1 << 3 ? 5 : 4 :
                            b < 1 << 6 ? b < 1 << 5 ? 3 : 2 :
                                         b < 1 << 7 ? 1 : 0;
    }
    private static int tz(byte b) // trailing zeros
    {
        return 7 - lz((byte)(b & -(sbyte)b));
    }
}

2012/09/28

Number of bits Fibonacci number

// Binet's formula: 
//
//         Fib(n) = ( Phi^n - (-1)^n / Phi^n ) / 5^(1/2)
//
// Where Phi = (1+5^(1/2))/2 ~ 1.61803398874989484820458683436563811772030917980576....(Knott)
//                           ~ 1.6180339887498948482045868343656381177203+             (Knuth)
//
// For larger values of n:
//
//         Fib(n) ~ Phi^n / 5^(1/2)
//
// For the number of bits, x, of Fib(n):
//
//            2^x ~ Fib(n) ~ Phi^n / 5^(1/2)
//
// Taking natural logarithms:
//
//        ln(2^x) ~ ln(Phi^n / 5^(1/2))
//
// So:
//      x * ln(2) ~ n * ln(Phi) - 1/2 * ln(5)
//              x ~ n * ln(Phi)/ln(2) - ln(5)/ln(2)/2
//
// Using "windows calculator" which seems to have a precision of ~105 bits ~31 decimal digits:
//
//              x ~ n * 0,6942419136306173017387902668986 - 1,1609640474436811739351597147447
//
// A partial bit (digit) is a bit (digit):
//
//        x ~ (int)(n * 0,6942419136306173017387902668986 - 1,1609640474436811739351597147447 + 1)
//
// The mantissa of a double has a precision of 52 bits ~ 16 decimal digits:
//
//        x = (int)(n * 0,6942419136306173 - 0,1609640474436812)
//
// Below it is correct for 0 <= n <= 24 * 10^6, after two days the output is: "R E A D Y"
// Also for n = 5 * 10^8 it's correct
//      for n = 10^9 it seems to be correct too.

using System;
using Xint = System.Numerics.BigInteger;
class bitLength_Fibonacci
{
    private static int bL_F(int n)
    {
        return n < 6 ? ++n / 2 : (int)(0.6942419136306173 * n - 0.1609640474436812);
    }

    static void Main()
    {
        Xint[] F = { 0, 1 };
        for (int n = 1; n <= 24000000; n++)
        {
            F = new Xint[2] { F[1], F[1] + F[0] };
            if (bL(F[0]) != bL_F(n))
            {
                Console.WriteLine(n + "   " + bL(F[0]) + "   " + bL_F(n));
            }
        }
        Console.WriteLine("R E A D Y");
        Console.ReadLine();
    }

    private static int bL(Xint X)
    {
        byte[] bytes = (X.Sign * X).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;
    }
}

2012/09/09

Fibonacci sequence binary plot, Ed Pegg Jr.(2003)

At Wolfram Mathworld you can find an intriguing plot of the
Fibonacci sequence represented in binary.
The one below shows it for small value's, 
Fibonacci(5)=5,  101 binary,        black, white, black    
Fibonacci(6)=8, 1000 binary, black, white, white, white

  




















A larger one: Fibonacci(1730)

Fibonacci(1000000000) Challenge in 0.1 ms

// Find the first 10 bytes and the last 10 bytes of f(1_000_000_000).
// That's right, fibonacci of one billion.

// I suppose the first 80 bits and the last 80 bits are more interesting,
// when you try to (dis)prove .... things.
// And, if you try, start with small numbers, first bit, last bit, first/last two bits, etc.
// By the way, the first bit is always one,
// though there is an exception, do you know it? (hint: start with small numbers)

// You can find the challenge at:
// weblogs.java.net/blog/kabutz/archive/2012/02/24/fibonacci-1000000000-challenge
// Or google "fibonacci 1000000000".

// My result(using 1 core):
//
//         F(1000000000)
//         bits: 694241913
//          MSB: 01 62 80 B8 2D 8C BE 0E DC 1B
//         time: 46 us
//          LSB: A9 53 2D F4 D2 D2 5B 5D B6 3B
//         time: 46 us

// Slightly less than 0.1 ms!!!!
// (the conversion of an 80 bits number to a string takes about 2 us, eternal question:
// can this conversion be done faster? I'm not going to answer it;)

// When the Most Significant Bytes are calculated with a precision of 80 bits,
// the lowest 2 or 3 bytes of these MSB's are wrong, if it's doubled to 160 bits,
// things look fine, are the same as those of Rex Young, which can be found at:
// weblogs.java.net/blog/kabutz/archive/2012/02/24/fibonacci-1000000000-challenge
// 
// Big question: "how to (dis)prove it?"
// It probably boils down to "more bits, higher accuracy, analyse the error". 
// Amazing: MSB's (multiplication of 160 bits numbers)
//      and LSB's (multiplication of  80 bits numbers)
//      take nearly the same time.
// I started with LowFib, returning the LSB's,
//   and thought HighFib, returning the MSB's,
//   would be 4 times slower.

// A few ways to speed it up:
// -upto a bitlength of 80 bits LowFib and HighFib do the same (iterations).
// -in the last iteration F(n) and F(n+1) are computed, get rid of F(n+1)
// -don't use a biginteger library, use arrays of uints,
//  and (Knuth's) classical multiplication algorithm.
// -translate it to c (or assembler)

// LSB's, 80 bits result from 80 bits multiplied by 80 bits.
// Looks like some time can be saved within the classical * algo.

// The final digits of fibonacci numbers repeat with a certain cycle length, see: 
// www.maths.surrey.ac.uk/hosted-sites/R.Knott/Fibonacci/fibmaths.html#cycles
// Looks like a modulo (remainder) operation.

// An underbound for the number of bits of fibonacci(n) is (2*n+1)/3
// Are there better bounds? 

// FibonacciHighLow(10^12)? The trillionth one.
// Use a few longs, it would surprise me if it takes more than 0.2 ms.

//       !!! Project >> Add Reference >> System.Numerics !!!
using Xint = System.Numerics.BigInteger;
using System.Diagnostics;
using System;

class FibonacciHighLow
{
    private static int nob;                           // number of bits
    private static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        for (int n = 1953125; n <= 1000000000; n *= 2)
        {
            Console.WriteLine("F(" + n + ")");
            string s = hexXint(HighFib(n));
            Console.WriteLine("bits: " + nob);
            Console.WriteLine(" MSB: " + s);
            HighFib(n);
            sw.Restart();
            for (int i = 0; i < 1000; i++)
            {
                HighFib(n);
                //hexXint(HighFib(n));               
            }
            sw.Stop();
            Console.WriteLine("time: " + sw.ElapsedMilliseconds + " us");
            s = hexXint(LowFib(n));
            Console.WriteLine(" LSB: " + s);
            LowFib(n);
            sw.Restart();
            for (int i = 0; i < 1000; i++)
            {
                LowFib(n);
                //hexXint(LowFib(n));
            }
            sw.Stop();
            Console.WriteLine("time: " + sw.ElapsedMilliseconds + " us");
            Console.WriteLine();
        }
        Console.ReadLine();
    }

    public static Xint LowFib(int n)
    {
        Xint Mask = (Xint.One << 80) - 1;
        Xint[] F = { 0, 1 };
        for (int i = fL2(n); i >= 0; i--)             // floorLog2
        {
            Xint L = 2 * F[1] - F[0];                 // L is a Lucas Number(2,1,3,4,7,..)
            switch ((3 << i & n) >> i)
            {
                case 0: F = new Xint[2] { L * F[0], L * F[1] - 1 }; break;
                case 2: F = new Xint[2] { L * F[0], L * F[1] + 1 }; break;
                case 1: F = new Xint[2] { L * F[1] - 1, L * (F[1] + F[0]) - 1 }; break;
                case 3: F = new Xint[2] { L * F[1] + 1, L * (F[1] + F[0]) + 1 }; break;
            }
            F[0] &= Mask;
            F[1] &= Mask;
        }
        return F[0];
    }
    // To do: 
    // case 0: F[1] = (L * F[1] - 1) & Mask;
    // if (F[1] < 0) F[1] = Mask; // ????
    // Looks like LowFib can be exact.
    // Test it with small Mask values (decimal 1,3,7,15,.../binary 1,11,111,1111,....)
    //    A few years ago I wrote Knuth's classical division algorithm in TWITS (TWo bITS)
    //    It's quite difficult to test it when bytes, ushorts or uints are used.
    //    With twits it's a lot easier.
    // Or, even better, don't, because:
    // mathworld.wolfram.com/FibonacciNumber.html
    // F(2^n + 2^(n+1)) ends in n+2 zeroes binary.
    // So 80 zeroes, n=78, 2^78+2^79, F(906.694.364.710.971.881.029.632) ends in 80 zeroes.

    public static Xint HighFib(int n)
    {
        Xint Mask = (Xint.One << 160) - 1;
        Xint[] F = { 0, 1 };
        int i = fL2(n);
        for (; i >= 0 && bL(F[0]) <= 160; i--)        // bitLength
        {
            Xint L = 2 * F[1] - F[0];
            switch ((3 << i & n) >> i)
            {
                case 0: F = new Xint[2] { L * F[0], L * F[1] - 1 }; break;
                case 2: F = new Xint[2] { L * F[0], L * F[1] + 1 }; break;
                case 1: F = new Xint[2] { L * F[1] - 1, L * (F[1] + F[0]) - 1 }; break;
                case 3: F = new Xint[2] { L * F[1] + 1, L * (F[1] + F[0]) + 1 }; break;
            }

        }
        int s = bL(F[0]) - 160;                       // shift
        int t = s;                                    // total shift
        F[0] >>= s;
        F[1] >>= s;
        for (; i >= 0; i--)
        {
            Xint L = 2 * F[1] - F[0];
            if ((1 << i & n) >> i == 0)
            {
                F[0] *= L; F[1] *= L;
            }
            else
            {
                F = new Xint[2] { L * F[1], L * (F[1] + F[0]) };
            }
            t *= 2;
            s = bL(F[0]) - 160;
            t += s;
            F[0] >>= s;
            F[1] >>= s;

        }
        t += bL(F[0]);
        nob = t;
        return F[0] >>= 80 + (8 - t & 7);
    }
    // mathworld.wolfram.com/FibonacciNumber.html
    // Edd Pegg Jr.'s plot(2003) might indicate that
    // repeating zeroes (binary) in Fibonacci Numbers,
    // in our case 160 zeroes, arise for F(VERY VERY LARGE)

    private static int fL2(int n) { int i = -1; while (n > 0) { i++; n /= 2; } return i; }

    private static int bL(Xint X)
    {
        byte[] bytes = X.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 string hexXint(Xint X)
    {
        byte[] bytes = X.ToByteArray();
        byte[] zeros = new byte[10];
        int i = bytes.Length;
        if (i == 11) i--;
        for (--i; i >= 0; i--)
        {
            zeros[i] = bytes[i];
        }
        string s = zeros[9].ToString("X2");
        for (i = 8; i >= 0; i--)
        {
            s += " " + zeros[i].ToString("X2");
        }
        return s;
    }
}

2012/02/24

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();
    }
}