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

/************************************************************
 * FUSS-CATALAN NUMBERS
 *
 * The nth Fuss-Catalan number is
 *
 *                 r       / pn + r \
 *   FC(n,p,r) = --------- (          ).
 *               pn + r    \    n    /
 *
 * For r = 1, the cases p = 2, 3, and 4 count
 * ordered full binary, ternary, and quaternary
 * trees with n internal nodes, respectively.
 *
 * We use the recurrence
 *
 *                     p(pn-p+1)...(pn-1)
 *   C_p(n) = C_p(n-1) ---------------------
 *                    ((p-1)n-p+3)...
 *                         ((p-1)n+1)
 *
 * with C_p(0) = 1.
 *
 * The numerator and denominator are factored.
 * We cancel common factors before multiplying,
 * preventing unnecessary intermediate overflow.
 *
 * We use uint64_t because unsigned long is not
 * guaranteed to be 64 bits on every platform.
 ************************************************************/

using number = std::uint64_t;

// Function declarations

int main();

number catalan2( unsigned int const n );
number catalan3( unsigned int const n );
number catalan4( unsigned int const n );

number fuss_catalan(
    unsigned int const n,
    unsigned int const p
);

number gcd(
    number m,
    number n
);

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

int main() {
    number const c2[]{
        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
    };

    number const c3[]{
        1ULL, 1ULL, 3ULL, 12ULL, 55ULL,
        273ULL, 1428ULL, 7752ULL, 43263ULL,
        246675ULL, 1430715ULL, 8414640ULL,
        50067108ULL, 300830572ULL, 1822766520ULL,
        11124755664ULL, 68328754959ULL,
        422030545335ULL, 2619631042665ULL,
        16332922290300ULL, 102240109897695ULL,
        642312451217745ULL, 4048514844039120ULL,
        25594403741131680ULL, 162250238001816900ULL,
        1031147983159782228ULL,
        6568517413771094628ULL
    };

    number const c4[]{
        1ULL, 1ULL, 4ULL, 22ULL, 140ULL,
        969ULL, 7084ULL, 53820ULL, 420732ULL,
        3362260ULL, 27343888ULL, 225568798ULL,
        1882933364ULL, 15875338990ULL,
        134993766600ULL, 1156393243320ULL,
        9969937491420ULL, 86445222719724ULL,
        753310723010608ULL, 6594154339031800ULL,
        57956002331347120ULL, 511238042454541545ULL,
        4524678117939182220ULL
    };

    std::size_t const C2{
        sizeof( c2 )/sizeof( c2[0] )
    };

    std::size_t const C3{
        sizeof( c3 )/sizeof( c3[0] )
    };

    std::size_t const C4{
        sizeof( c4 )/sizeof( c4[0] )
    };

    std::cout << "Binary Catalan numbers:"
              << std::endl;

    for (
        unsigned int n{ 0 };
        n < C2;
        ++n
    ) {
        number const result{ catalan2( n ) };

        assert( result == c2[n] );

        std::cout << "C2(" << n << ") = "
                  << result << std::endl;
    }

    std::cout << std::endl
              << "Ternary Catalan numbers:"
              << std::endl;

    for (
        unsigned int n{ 0 };
        n < C3;
        ++n
    ) {
        number const result{ catalan3( n ) };

        assert( result == c3[n] );

        std::cout << "C3(" << n << ") = "
                  << result << std::endl;
    }

    std::cout << std::endl
              << "Quaternary Catalan numbers:"
              << std::endl;

    for (
        unsigned int n{ 0 };
        n < C4;
        ++n
    ) {
        number const result{ catalan4( n ) };

        assert( result == c4[n] );

        std::cout << "C4(" << n << ") = "
                  << result << std::endl;
    }

    std::cout << std::endl
              << "All tests passed."
              << std::endl;

    return 0;
}

/************************************************************
 * BINARY CATALAN NUMBERS
 *
 *                 2(2n - 1)
 *   C2(n) = ------------------- C2(n - 1)
 *                    n + 1
 *
 * Supported for 0 <= n <= 36.
 ************************************************************/

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

    return fuss_catalan( n, 2 );
}

/************************************************************
 * TERNARY CATALAN NUMBERS
 *
 *              3(3n - 2)(3n - 1)
 *   C3(n) = ----------------------- C3(n - 1)
 *                (2n)(2n + 1)
 *
 * Supported for 0 <= n <= 26.
 ************************************************************/

number catalan3(
    unsigned int const n
) {
    assert( n <= 26 );

    return fuss_catalan( n, 3 );
}

/************************************************************
 * QUATERNARY CATALAN NUMBERS
 *
 *          4(4n - 3)(4n - 2)(4n - 1)
 * C4(n) = ---------------------------- C4(n - 1)
 *          (3n - 1)(3n)(3n + 1)
 *
 * Supported for 0 <= n <= 22.
 ************************************************************/

number catalan4(
    unsigned int const n
) {
    assert( n <= 22 );

    return fuss_catalan( n, 4 );
}

/************************************************************
 * GENERAL FUSS-CATALAN CALCULATION FOR r = 1
 *
 * Calculate successive values C_p(1), ..., C_p(n)
 * using the factored recurrence.
 *
 * Each recurrence step has p numerator factors
 * and p - 1 denominator factors.
 *
 * For p = 3, these are:
 *
 *   Numerator:   3, 3k - 2, 3k - 1
 *
 *   Denominator: 2k, 2k + 1
 *
 * For p = 4:
 *
 *   Numerator:   4, 4k - 3, 4k - 2, 4k - 1
 *
 *   Denominator: 3k - 1, 3k, 3k + 1
 *
 * For each denominator factor:
 *
 *   1. Cancel common factors with every
 *      numerator factor.
 *
 *   2. Cancel any remaining factors with
 *      the previous Catalan number.
 *
 * Once all denominator factors are cancelled,
 * multiply the remaining numerator factors.
 *
 * This order avoids multiplying by a large
 * numerator before performing the division.
 *
 * Preconditions:
 *   2 <= p <= 4
 *   n is within the supported range for p.
 ************************************************************/

number fuss_catalan(
    unsigned int const n,
    unsigned int const p
) {
    assert( p >= 2 );
    assert( p <= 4 );

    assert(
        ( p == 2 && n <= 36 )
        || ( p == 3 && n <= 26 )
        || ( p == 4 && n <= 22 )
    );

    number result{ 1 };

    for (
        unsigned int k{ 1 };
        k <= n;
        ++k
    ) {
        // Four numerator factors and three
        // denominator factors are sufficient
        // for all three supported cases.
        number numerator[4]{};
        number denominator[3]{};

        numerator[0] = p;

        for (
            unsigned int j{ 1 };
            j < p;
            ++j
        ) {
            numerator[j] =
                static_cast<number>( p )*k - p + j;
        }

        for (
            unsigned int j{ 0 };
            j < p - 1;
            ++j
        ) {
            denominator[j] =
                static_cast<number>( p - 1 )*k
                - p + 3 + j;
        }

        // Cancel every denominator factor.
        for (
            unsigned int j{ 0 };
            j < p - 1;
            ++j
        ) {
            // First cancel against the
            // numerator factors.
            for (
                unsigned int i{ 0 };
                i < p;
                ++i
            ) {
                number const divisor{
                    gcd(
                        numerator[i],
                        denominator[j]
                    )
                };

                numerator[i] /= divisor;
                denominator[j] /= divisor;
            }

            // Cancel any remaining denominator
            // against the preceding result.
            number const divisor{
                gcd( result, denominator[j] )
            };

            result /= divisor;
            denominator[j] /= divisor;

            // The recurrence must produce
            // an integer.
            assert( denominator[j] == 1 );
        }

        // All denominators have been cancelled.
        // Multiply the remaining numerator factors.
        for (
            unsigned int i{ 0 };
            i < p;
            ++i
        ) {
            assert( numerator[i] != 0 );

            assert(
                result
                <= std::numeric_limits<number>::max()
                   /numerator[i]
            );

            result *= numerator[i];
        }
    }

    return result;
}

/************************************************************
 * GREATEST COMMON DIVISOR
 *
 * Calculate gcd(m, n) using the Euclidean algorithm.
 *
 * The result is positive unless both arguments
 * are zero.
 *
 * In this program, gcd() is called only with
 * positive factors, so its result is always
 * at least one.
 ************************************************************/

number gcd(
    number m,
    number n
) {
    assert( m != 0 || n != 0 );

    while ( n != 0 ) {
        number const remainder{ m%n };

        m = n;
        n = remainder;
    }

    return m;
}