#include <iostream>
#include <cstdint>
#include <cstddef>
#include <limits>
#include <numeric>
#include <cassert>

/************************************************************
 * CATALAN NUMBERS
 *
 * The nth Catalan number counts the paths on an n-by-n
 * square grid from the top-left corner to the bottom-right
 * corner such that:
 *
 *   1. Each move is one step right or one step down.
 *
 *   2. The path never crosses below the main diagonal.
 *
 * Equivalently, at every point along the path, the number
 * of rightward moves is at least the number of downward
 * moves.
 *
 * For n = 3, there are five such paths:
 *
 *     RRRDDD
 *     RRDRDD
 *     RRDDRD
 *     RDRRDD
 *     RDRDRD
 *
 * Thus C(3) = 5.
 *
 * We will calculate Catalan numbers using several methods:
 *
 *   1. A binomial coefficient divided by n + 1.
 *   2. The difference of two binomial coefficients.
 *   3. A formula involving three factorials.
 *   4. A product of rational factors.
 *   5. The sum recurrence using naive recursion.
 *   6. The sum recurrence using dynamic programming.
 *   7. The ratio recurrence using recursion.
 *   8. The ratio recurrence using iteration.
 *
 * Some formulas overflow even when the final Catalan
 * number fits in a 64-bit unsigned integer.
 *
 * Consequently, the demonstration restricts each
 * implementation to a range in which its intermediate
 * calculations are safe.
 ************************************************************/

using number = std::uint64_t;

// Function declarations

int main();

number catalan_bin( unsigned int const n );
number catalan_dif( unsigned int const n );
number catalan_fac( unsigned int const n );
number catalan_prd( unsigned int const n );
number catalan_sum( unsigned int const n );
number catalan_dyn( unsigned int const n );
number catalan_rec( unsigned int const n );
number catalan_opt( unsigned int const n );

number factorial( unsigned int const n );

number binomial(
    unsigned int const n,
    unsigned int const k
);

// Calculate C(k) from C(k - 1), cancelling
// factors before multiplying.
number catalan_next(
    number       const previous,
    unsigned int const k
);

/************************************************************
 * MAIN
 ************************************************************/

int main() {
    number const expected[]{
        1ULL,
        1ULL,
        2ULL,
        5ULL,
        14ULL,
        42ULL,
        132ULL,
        429ULL,
        1430ULL,
        4862ULL,
        16796ULL,
        58786ULL,
        208012ULL,
        742900ULL,
        2674440ULL,
        9694845ULL,
        35357670ULL,
        129644790ULL,
        477638700ULL,
        1767263190ULL,
        6564120420ULL,
        24466267020ULL,
        91482563640ULL,
        343059613650ULL,
        1289904147324ULL,
        4861946401452ULL,
        18367353072152ULL,
        69533550916004ULL,
        263747951750360ULL,
        1002242216651368ULL,
        3814986502092304ULL,
        14544636039226909ULL,
        55534064877048198ULL,
        212336130412243110ULL,
        812944042149730764ULL,
        3116285494907301262ULL,
        11959798385860453492ULL
    };

    std::size_t const count{
        sizeof( expected )/sizeof( expected[0] )
    };

    for (
        unsigned int n{ 0 };
        n < count;
        ++n
    ) {
        std::cout << "C(" << n << ") = "
                  << expected[n] << std::endl;

        // The central binomial coefficient
        // fits in 64 bits through n = 33.
        if ( n <= 33 ) {
            number const value{ catalan_bin( n ) };

            assert( value == expected[n] );

            std::cout << "Binomial divided by n + 1: "
                      << value << std::endl;

            number const difference{
                catalan_dif( n )
            };

            assert( difference == expected[n] );

            std::cout << "Difference of binomials:   "
                      << difference << std::endl;
        }

        // The factorial and direct-product
        // demonstrations use smaller ranges
        // to avoid intermediate overflow.
        if ( n <= 10 ) {
            number const factorial_value{
                catalan_fac( n )
            };

            number const product_value{
                catalan_prd( n )
            };

            assert( factorial_value == expected[n] );
            assert( product_value == expected[n] );

            std::cout << "Factorial formula:         "
                      << factorial_value << std::endl;

            std::cout << "Product formula:           "
                      << product_value << std::endl;
        }

        // Naive recursion becomes expensive
        // surprisingly quickly.
        if ( n <= 12 ) {
            number const value{ catalan_sum( n ) };

            assert( value == expected[n] );

            std::cout << "Naive recursive sum:       "
                      << value << std::endl;
        }

        // These implementations can calculate
        // every value in the demonstration.
        number const dynamic_value{
            catalan_dyn( n )
        };

        number const recursive_value{
            catalan_rec( n )
        };

        number const iterative_value{
            catalan_opt( n )
        };

        assert( dynamic_value == expected[n] );
        assert( recursive_value == expected[n] );
        assert( iterative_value == expected[n] );

        std::cout << "Dynamic programming:       "
                  << dynamic_value << std::endl;

        std::cout << "Recursive ratio:           "
                  << recursive_value << std::endl;

        std::cout << "Iterative ratio:           "
                  << iterative_value << std::endl;

        std::cout << std::endl;
    }

    return 0;
}

/************************************************************
 * METHOD 1: BINOMIAL COEFFICIENT
 *
 *                1       (2n)!
 *     C(n) = --------- * -------
 *              n + 1     n! n!
 *
 * Equivalently,
 *
 *     C(n) = binomial(2n, n)/(n + 1).
 *
 * The intermediate binomial coefficient must
 * fit in the chosen integer type.
 *
 * Supported here for 0 <= n <= 33.
 ************************************************************/

number catalan_bin(
    unsigned int const n
) {
    assert( n <= 33 );

    return binomial( 2*n, n )/(n + 1);
}

/************************************************************
 * METHOD 2: DIFFERENCE OF BINOMIAL COEFFICIENTS
 *
 *     C(n) = binomial(2n, n)
 *          - binomial(2n, n + 1).
 *
 * Both binomial coefficients must fit in
 * the chosen integer type.
 *
 * Supported here for 0 <= n <= 33.
 ************************************************************/

number catalan_dif(
    unsigned int const n
) {
    assert( n <= 33 );

    return binomial( 2*n, n )
         - binomial( 2*n, n + 1 );
}

/************************************************************
 * METHOD 3: FACTORIAL FORMULA
 *
 *                (2n)!
 *     C(n) = -------------
 *             n! (n + 1)!
 *
 * Factorials grow very quickly.
 *
 * Since 20! fits in a 64-bit unsigned integer
 * but 21! does not, this implementation is
 * restricted to n <= 10.
 ************************************************************/

number catalan_fac(
    unsigned int const n
) {
    assert( n <= 10 );

    return factorial( 2*n )
         / factorial( n )
         / factorial( n + 1 );
}

/************************************************************
 * METHOD 4: PRODUCT FORMULA
 *
 *             n
 *             ---
 *             | |  (n + k)/k
 *             ---
 *            k = 2
 *
 * We multiply the numerators and denominators
 * separately, then perform one division.
 *
 * This avoids truncating individual rational
 * factors, but the products grow rapidly.
 *
 * Supported here for 0 <= n <= 10.
 ************************************************************/

number catalan_prd(
    unsigned int const n
) {
    assert( n <= 10 );

    number numerator{ 1 };
    number denominator{ 1 };

    for (
        unsigned int k{ 2 };
        k <= n;
        ++k
    ) {
        numerator *= n + k;
        denominator *= k;
    }

    assert( numerator%denominator == 0 );

    return numerator/denominator;
}

/************************************************************
 * METHOD 5: NAIVE RECURSION
 *
 *                   n - 1
 *                    ---
 *                    \
 *     C(n) =          /   C(k)*C(n - k - 1)
 *                    ---
 *                   k = 0
 *
 * with C(0) = 1.
 *
 * This formula follows naturally from splitting
 * a Catalan structure into two smaller structures.
 *
 * Unfortunately, naive recursion calculates
 * the same Catalan numbers repeatedly.
 *
 * This implementation is restricted to n <= 12
 * in the demonstration to keep execution quick.
 ************************************************************/

number catalan_sum(
    unsigned int const n
) {
    assert( n <= 12 );

    if ( n == 0 ) {
        return 1;
    }

    number result{ 0 };

    for (
        unsigned int k{ 0 };
        k < n;
        ++k
    ) {
        result += catalan_sum( k )
                * catalan_sum( n - k - 1 );
    }

    return result;
}

/************************************************************
 * METHOD 6: DYNAMIC PROGRAMMING
 *
 * Use the same recurrence as catalan_sum(),
 * but calculate each Catalan number only once.
 *
 * Start with C[0] = 1.
 *
 * Then calculate C[1], C[2], ..., C[n],
 * storing each result for later use.
 *
 * Unlike a variable-length array, a dynamically
 * allocated array is valid standard C++.
 *
 * The caller of new[] must eventually use delete[].
 *
 * Supported here for 0 <= n <= 36.
 ************************************************************/

number catalan_dyn(
    unsigned int const n
) {
    assert( n <= 36 );

    number *C{ new number[n + 1]{} };

    C[0] = 1;

    for (
        unsigned int current{ 1 };
        current <= n;
        ++current
    ) {
        for (
            unsigned int k{ 0 };
            k < current;
            ++k
        ) {
            C[current] +=
                C[k]*C[current - k - 1];
        }
    }

    number const result{ C[n] };

    delete[] C;

    return result;
}

/************************************************************
 * HELPER: RATIO RECURRENCE
 *
 *                 2(2k - 1)
 *     C(k) = ------------------- C(k - 1)
 *                   k + 1
 *
 * A direct calculation might overflow before
 * the division, even if C(k) fits.
 *
 * Instead, cancel common factors between the
 * numerator, denominator, and previous result
 * before multiplying.
 *
 * The recurrence is exact, so the denominator
 * must eventually reduce to one.
 *
 * Supported here for 1 <= k <= 36.
 ************************************************************/

number catalan_next(
    number       const previous,
    unsigned int const k
) {
    assert( k >= 1 );
    assert( k <= 36 );

    number result{ previous };

    number numerator{ 2*(2*k - 1) };
    number denominator{ k + 1 };

    // First cancel factors from the
    // previous Catalan number.
    number divisor{
        std::gcd( result, denominator )
    };

    result /= divisor;
    denominator /= divisor;

    // Cancel any remaining denominator
    // against the numerator.
    divisor = std::gcd( numerator, denominator );

    numerator /= divisor;
    denominator /= divisor;

    // The recurrence produces an integer.
    assert( denominator == 1 );

    // Check that the final multiplication
    // fits in our integer type.
    assert(
        result
        <= std::numeric_limits<number>::max()/numerator
    );

    return result*numerator;
}

/************************************************************
 * METHOD 7: RECURSIVE RATIO
 *
 * Use the ratio recurrence recursively.
 *
 * Each call calculates C(n - 1) once.
 *
 * Unlike the naive sum recurrence, this
 * implementation makes only n recursive calls.
 *
 * Supported here for 0 <= n <= 36.
 ************************************************************/

number catalan_rec(
    unsigned int const n
) {
    assert( n <= 36 );

    if ( n == 0 ) {
        return 1;
    }

    return catalan_next(
        catalan_rec( n - 1 ),
        n
    );
}

/************************************************************
 * METHOD 8: ITERATIVE RATIO
 *
 * Use the same ratio recurrence, but calculate
 * successive Catalan numbers in a loop.
 *
 * We need to store only the current result.
 *
 * Supported here for 0 <= n <= 36.
 ************************************************************/

number catalan_opt(
    unsigned int const n
) {
    assert( n <= 36 );

    number result{ 1 };

    for (
        unsigned int k{ 1 };
        k <= n;
        ++k
    ) {
        result = catalan_next( result, k );
    }

    return result;
}

/************************************************************
 * FACTORIAL
 *
 * Calculate n!.
 *
 * Since 20! is the largest factorial that
 * fits in a 64-bit unsigned integer,
 * require n <= 20.
 ************************************************************/

number factorial(
    unsigned int const n
) {
    assert( n <= 20 );

    number result{ 1 };

    for (
        unsigned int k{ 2 };
        k <= n;
        ++k
    ) {
        result *= k;
    }

    return result;
}

/************************************************************
 * BINOMIAL COEFFICIENT
 *
 * Calculate binomial(n, k) without calculating
 * three factorials.
 *
 * Use symmetry to replace k with min(k, n - k).
 *
 * Then calculate the product
 *
 *     n - k + 1     n - k + 2             n
 *     ---------  *  ---------  * ... *  ---
 *         1             2                 k
 *
 * At each step, cancel common factors before
 * multiplying. This reduces the risk of
 * intermediate overflow.
 *
 * If k > n, return zero.
 *
 * The caller must still ensure that the
 * binomial coefficient fits in 64 bits.
 ************************************************************/

number binomial(
    unsigned int const n,
    unsigned int const k
) {
    assert( n <= 66 );

    if ( k > n ) {
        return 0;
    }

    unsigned int const smaller{
        ( k < n - k ) ? k : n - k
    };

    number result{ 1 };

    for (
        unsigned int i{ 1 };
        i <= smaller;
        ++i
    ) {
        number numerator{ n - smaller + i };
        number denominator{ i };

        // Cancel factors within the
        // current numerator and denominator.
        number divisor{
            std::gcd( numerator, denominator )
        };

        numerator /= divisor;
        denominator /= divisor;

        // Cancel the remaining denominator
        // against the accumulated result.
        divisor = std::gcd( result, denominator );

        result /= divisor;
        denominator /= divisor;

        // Each intermediate binomial
        // coefficient is an integer.
        assert( denominator == 1 );

        // Avoid overflowing the multiplication.
        assert(
            result
            <= std::numeric_limits<number>::max()/numerator
        );

        result *= numerator;
    }

    return result;
}