#include <iostream>
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <stdexcept>
#include <vector>

// Seven methods for calculating binomial coefficients.
// Requires C++11 or later. All arithmetic uses 64-bit unsigned integers.
using number = std::uint64_t;
number const LIMIT{ std::numeric_limits<number>::max() };

number factorial( unsigned int const n );
number gcd( number m, number n );
number checked_add( number const x, number const y );
number checked_multiply( number const x, number const y );
number binomial_fac( unsigned int const n, unsigned int const k );
number binomial_pr1( unsigned int const n, unsigned int const k );
number binomial_pr2( unsigned int const n, unsigned int const k );
number binomial_gcd( unsigned int const n, unsigned int const k );
number binomial_rec( unsigned int const n, unsigned int const k );
number binomial_dy1( unsigned int const n, unsigned int const k );
number binomial_dy2( unsigned int const n, unsigned int const k );

// Print a result or explain why this particular method cannot calculate it.
void demonstrate( char const *label, number (*method)(unsigned int, unsigned int),
                  unsigned int const n, unsigned int const k ) {
    std::cout << label << ": ";
    try {
        std::cout << method( n, k );
    } catch ( std::range_error const &error ) {
        std::cout << "unavailable (" << error.what() << ")";
    }
    std::cout << std::endl;
}

int main() {
    number const binomial_100_k[]{
        1ULL, 100ULL, 4950ULL, 161700ULL, 3921225ULL,
        75287520ULL, 1192052400ULL, 16007560800ULL,
        186087894300ULL, 1902231808400ULL, 17310309456440ULL,
        141629804643600ULL, 1050421051106700ULL,
        7110542499799200ULL, 44186942677323600ULL,
        253338471349988640ULL, 1345860629046814650ULL,
        6650134872937201800ULL
    };
    number const binomial_281_k[]{
        1ULL, 281ULL, 39340ULL, 3658620ULL, 254274090ULL,
        14086784586ULL, 647992090956ULL, 25456832144700ULL,
        871896500955975ULL, 26447527195664575ULL,
        719372739722076440ULL, 17722728405880246840ULL
    };

    assert( gcd( 0, 0 ) == 0 );
    assert( gcd( 12, 18 ) == 6 );
    assert( gcd( 0, 12 ) == 12 );

    // All seven methods must handle the boundary cases.
    number (*methods[])(unsigned int, unsigned int){
        binomial_fac, binomial_pr1, binomial_pr2, binomial_gcd,
        binomial_rec, binomial_dy1, binomial_dy2
    };
    for ( auto method : methods ) {
        assert( method( 0, 0 ) == 1 );
        assert( method( 5, 0 ) == 1 );
        assert( method( 5, 5 ) == 1 );
        assert( method( 5, 6 ) == 0 );
    }

    std::cout << "Binomial coefficients for n = 42" << std::endl;
    for ( unsigned int k{ 0 }; k <= 43; ++k ) {
        std::cout << "(42, " << k << ")" << std::endl;
        demonstrate( "  Factorial", binomial_fac, 42, k );
        demonstrate( "  First product", binomial_pr1, 42, k );
        demonstrate( "  Symmetric product", binomial_pr2, 42, k );
        demonstrate( "  GCD product", binomial_gcd, 42, k );
        if ( k > 42 || std::min( k, 42U - k ) <= 4 ) {
            demonstrate( "  Naive recursion", binomial_rec, 42, k );
        } else {
            std::cout << "  Naive recursion: skipped (too slow)" << std::endl;
        }
        demonstrate( "  Full Pascal triangle", binomial_dy1, 42, k );
        demonstrate( "  Reduced Pascal triangle", binomial_dy2, 42, k );
    }

    std::cout << "\nBinomial coefficients for n = 100" << std::endl;
    for ( unsigned int k{ 0 }; k < 18; ++k ) {
        number const exact{ binomial_100_k[k] };
        assert( binomial_gcd( 100, k ) == exact );
        assert( binomial_dy2( 100, k ) == exact );
        std::cout << "(100, " << k << ") = " << exact << std::endl;
        demonstrate( "  Factorial", binomial_fac, 100, k );
        demonstrate( "  First product", binomial_pr1, 100, k );
        demonstrate( "  Symmetric product", binomial_pr2, 100, k );
        demonstrate( "  GCD product", binomial_gcd, 100, k );
        if ( k <= 3 ) {
            demonstrate( "  Naive recursion", binomial_rec, 100, k );
        }
        demonstrate( "  Full Pascal triangle", binomial_dy1, 100, k );
        demonstrate( "  Reduced Pascal triangle", binomial_dy2, 100, k );
    }

    std::cout << "\nBinomial coefficients for n = 281" << std::endl;
    for ( unsigned int k{ 0 }; k < 12; ++k ) {
        number const exact{ binomial_281_k[k] };
        assert( binomial_gcd( 281, k ) == exact );
        assert( binomial_dy2( 281, k ) == exact );
        std::cout << "(281, " << k << ") = " << exact << std::endl;
    }

    // A central coefficient that fits and one that does not.
    assert( binomial_gcd( 67, 33 ) == 14226520737620288370ULL );
    assert( binomial_dy1( 67, 33 ) == 14226520737620288370ULL );
    assert( binomial_dy2( 67, 33 ) == 14226520737620288370ULL );
    bool overflow_detected{ false };
    try {
        (void)binomial_gcd( 68, 34 );
    } catch ( std::range_error const & ) {
        overflow_detected = true;
    }
    assert( overflow_detected );
    std::cout << "\nAll assertions passed." << std::endl;
}

number checked_add( number const x, number const y ) {
    if ( x > LIMIT - y ) {
        throw std::range_error{ "64-bit addition overflow" };
    }
    return x + y;
}

number checked_multiply( number const x, number const y ) {
    if ( x != 0 && y > LIMIT/x ) {
        throw std::range_error{ "64-bit multiplication overflow" };
    }
    return x*y;
}

number factorial( unsigned int const n ) {
    number result{ 1 };
    for ( unsigned int i{ 2 }; i <= n; ++i ) {
        result = checked_multiply( result, i );
    }
    return result;
}

// Method 1: n!/(k!(n-k)!). It may overflow before division.
number binomial_fac( unsigned int const n, unsigned int const k ) {
    if ( k > n ) return 0;
    number const numerator{ factorial( n ) };
    number const denominator{
        checked_multiply( factorial( k ), factorial( n - k ) )
    };
    return numerator/denominator;
}

// Method 2: multiply k numerator and denominator factors separately.
number binomial_pr1( unsigned int const n, unsigned int const k ) {
    if ( k > n ) return 0;
    number numerator{ 1 };
    number denominator{ 1 };
    for ( unsigned int i{ 1 }; i <= k; ++i ) {
        numerator = checked_multiply( numerator, n - i + 1 );
        denominator = checked_multiply( denominator, i );
    }
    return numerator/denominator;
}

// Method 3: use symmetry to reduce the number of factors.
number binomial_pr2( unsigned int const n, unsigned int const k ) {
    if ( k > n ) return 0;
    unsigned int const smaller{ std::min( k, n - k ) };
    return binomial_pr1( n, smaller );
}

// Euclid's algorithm. By convention, gcd(0,0) = 0.
number gcd( number m, number n ) {
    while ( n != 0 ) {
        number const remainder{ m%n };
        m = n;
        n = remainder;
    }
    return m;
}

// Method 4: exact cancellation at each step.
// At iteration i, result is C(n-k+i-1, i-1).
// C(n-k+i, i) = result*(n-k+i)/i, which is an integer.
// Cancelling i against result and the new numerator before
// multiplying ensures the intermediate product equals the result.
number binomial_gcd( unsigned int const n, unsigned int const k ) {
    if ( k > n ) return 0;
    unsigned int const smaller{ std::min( k, n - k ) };
    number result{ 1 };
    for ( unsigned int i{ 1 }; i <= smaller; ++i ) {
        number numerator{ static_cast<number>( n ) - smaller + i };
        number denominator{ i };
        number divisor{ gcd( result, denominator ) };
        result /= divisor;
        denominator /= divisor;
        divisor = gcd( numerator, denominator );
        numerator /= divisor;
        denominator /= divisor;
        assert( denominator == 1 );
        result = checked_multiply( result, numerator );
    }
    return result;
}

// Method 5: Pascal's identity without memoization; exponential time.
number binomial_rec( unsigned int const n, unsigned int const k ) {
    if ( k > n ) return 0;
    if ( k == 0 || k == n ) return 1;
    return checked_add(
        binomial_rec( n - 1, k ),
        binomial_rec( n - 1, k - 1 )
    );
}

// Method 6: build the entire row of Pascal's triangle in place.
// If ANY entry overflows, this full-row method cannot continue,
// even if the particular requested coefficient fits.
number binomial_dy1( unsigned int const n, unsigned int const k ) {
    if ( k > n ) return 0;
    if ( n > 67 ) {
        throw std::range_error{ "full Pascal row may overflow" };
    }
    std::vector<number> row( static_cast<std::size_t>( n ) + 1, 0 );
    row[0] = 1;
    for ( unsigned int m{ 1 }; m <= n; ++m ) {
        row[m] = 1;
        for ( unsigned int i{ m - 1 }; i > 0; --i ) {
            row[i] = checked_add( row[i], row[i - 1] );
        }
    }
    return row[k];
}

// Method 7: retain only coefficients up to min(k,n-k).
// Uses O(k) storage and O(n*k) additions; checks every addition.
number binomial_dy2( unsigned int const n, unsigned int const k ) {
    if ( k > n ) return 0;
    unsigned int const smaller{ std::min( k, n - k ) };
    std::vector<number> row( static_cast<std::size_t>( smaller ) + 1, 0 );
    row[0] = 1;
    for ( unsigned int m{ 1 }; m <= n; ++m ) {
        unsigned int const last{ std::min( m, smaller ) };
        for ( unsigned int i{ last }; i > 0; --i ) {
            row[i] = checked_add( row[i], row[i - 1] );
        }
    }
    return row[smaller];
}