python nCr

import math
import timeit
import typing
from fractions import Fraction
from functools import reduce
from math import factorial
from operator import mul  # or mul=lambda x,y:x*y

import gmpy2
import scipy
import sympy
from scipy.special import binom


def nasbanov(n, k):
    # https://stackoverflow.com/a/3027128/9044183
    return int(reduce(mul, (Fraction(n - i, i + 1) for i in range(k)), 1))


def andrewdalke(n, k):
    # https://stackoverflow.com/a/3025547/9044183
    if 0 <= k <= n:
        ntok = 1
        ktok = 1
        for t in range(1, min(k, n - k) + 1):
            ntok *= n
            ktok *= t
            n -= 1
        return ntok // ktok
    else:
        return 0


def jfsmodifiedscipy(N, k):
    # https://stackoverflow.com/a/3025547/9044183
    if (k > N) or (N < 0) or (k < 0):
        return 0
    top = N
    val = 1
    while top > (N - k):
        val *= top
        top -= 1
    n = 1
    while n < k + 1:
        val /= n
        n += 1
    return val


def stdfactorial(n, r):
    # https://stackoverflow.com/a/19109792/9044183
    return factorial(n) // factorial(r) // factorial(n - r)


def stdfactorial2(n, r):
    # https://stackoverflow.com/a/52987317/9044183
    return factorial(n) // (factorial(r) * factorial(n - r))


def recursive(n, k):
    # https://stackoverflow.com/a/43107785/9044183
    return 1 if k == 0 or k == n else recursive(n - 1, k - 1) + recursive(n - 1, k)


def pantelis300(n, k):
    # https://stackoverflow.com/a/25813723/9044183
    c = [0] * (n + 1)
    c[0] = 1
    for i in range(1, n + 1):
        c[i] = 1
        j = i - 1
        while j > 0:
            c[j] += c[j - 1]
            j -= 1

    return c[k]


def rojas(n, k):
    # https://stackoverflow.com/a/56120255/9044183
    m = 0
    if k == 0:
        m = 1
    if k == 1:
        m = n
    if k >= 2:
        num, dem, op1, op2 = 1, 1, k, n
        while op1 >= 1:
            num *= op2
            dem *= op1
            op1 -= 1
            op2 -= 1
        m = num // dem
    return m


def wirawan(n, r):
    # https://stackoverflow.com/a/12115521/9044183
    assert n >= 0
    assert 0 <= r <= n

    c = 1
    for (num, denom) in zip(range(n, n - r, -1), range(1, r + 1, 1)):
        c = (c * num) // denom
    return c


def kta(n, r):
    # https://stackoverflow.com/a/58247588/9044183
    if r == 0:
        return 1
    p = n
    for i in range(1, r):
        p = p * (n - i) / i
    else:
        p = p / r
    return p


def rabih(n, k):
    # https://stackoverflow.com/a/30428134/9044183
    if k == n:
        return 1
    if k > n:
        return 0
    d, q = max(k, n - k), min(k, n - k)
    num = 1
    for n in range(d + 1, n + 1):
        num *= n
    denom = 1
    for d in range(1, q + 1):
        denom *= d
    return num / denom


def scipybinom(n, r):
    return scipy.special.binom(n, r)


def scipycombexact(n, r):
    return scipy.special.comb(n, r, True)


def hybrid(n: typing.Union[int, float], k: typing.Union[int, float]) -> typing.Union[int, float]:
    # my own custom hybrid solution
    def is_integer(n):
        return isinstance(n, int) or n.is_integer()
    if k < 0:
        raise ValueError("k cannot be negative.")
    elif n == 0:
        return 0
    elif k == 0 or k == n:
        return 1
    elif is_integer(n) and is_integer(k):
        return int(gmpy2.comb(int(n), int(k)))
    elif n > 0:
        return scipy.special.binom(n, k)
    else:
        return float(sympy.binomial(n, k))


def test_basic_func(func):
    [[func(i, j) for j in range(i + 1)] for i in range(1, 100)]


def oobf(func):
    [func(100, j) for j in range(100, 200)]


def ffrac_n(func):
    [[func(i + 0.5, j) for j in range(i + 1)] for i in range(1, 100)]


def ffrac_r(func):
    [[func(i, j + 0.5) for j in range(i + 1)] for i in range(1, 100)]


def ffrac_nr(func):
    [[func(i + 0.5, j + 0.5) for j in range(i + 1)] for i in range(1, 100)]


def negative_test(func):
    [[func(-i, j) for j in range(i + 1)] for i in range(1, 100)]


def negative_fracnr_test(func):
    [[func(-i - 0.5, j + 0.5) for j in range(i + 1)] for i in range(1, 100)]


def test():
    # pos int n/r and negative n and out of bounds and frac n and frac r
    fully_featured = [sympy.binomial, hybrid]
    negative_n = [nasbanov, rojas, kta, gmpy2.comb] + fully_featured
    frac_n_and_r = [scipybinom] + fully_featured
    frac_n = [kta] + frac_n_and_r
    frac_r = [scipy.special.comb] + frac_n_and_r
    # large indices return 0
    out_of_bounds = [nasbanov, andrewdalke, jfsmodifiedscipy, rojas, kta, rabih, math.comb, gmpy2.comb, scipybinom,
                     scipy.special.comb, sympy.binomial, scipycombexact, hybrid]
    basicfuncs = [stdfactorial, stdfactorial2, pantelis300, wirawan] + out_of_bounds
    tests = {
        "basic": (basicfuncs, test_basic_func,),
        "oob": (out_of_bounds, oobf,),
        "fn": (frac_n, ffrac_n,),
        "fr": (frac_r, ffrac_r,),
        "fnr": (frac_n_and_r, ffrac_nr,),
        "nn": (negative_n, negative_test,),
        "f": (fully_featured, negative_fracnr_test),
    }
    for name, (funcs, testfunc) in tests.items():
        print(name.upper())
        functimes = {}
        for func in funcs:
            funcname = f"{getattr(func, '__module__', '?')}.{func.__name__}"
            print(f"timing {funcname}...")
            functime = timeit.timeit(lambda: testfunc(func), number=10) / 10
            functimes[funcname] = functime
        print("-----")
        functimes = {k: v for k, v in sorted(functimes.items(), key=lambda item: item[1])}
        for k, v in functimes.items():
            print(k, v)


def get_features():
    funcs = [nasbanov, andrewdalke, jfsmodifiedscipy, stdfactorial, stdfactorial2, recursive,
             pantelis300, rojas, wirawan, kta, rabih, math.comb, gmpy2.comb, scipybinom,
             scipy.special.comb, sympy.binomial, scipycombexact]
    for func in funcs:
        funcname = f"{getattr(func, '__module__', '?')}.{func.__name__}"
        features = ["pos int n/r"]
        try:
            assert func(-4, 2) == 10
            features.append("negative n")
        except Exception as e:
            pass
        try:
            assert func(5, 10) == 0
            features.append("out of bounds")
        except Exception as e:
            pass
        try:
            assert func(0.5, 3) == 1 / 16
            features.append("frac n")
        except Exception as e:
            pass
        try:
            assert not float(func(3, 0.5)).is_integer()
            features.append("frac r")
        except Exception as e:
            pass

        print(f"{funcname} supports {' and '.join(features)}")


if __name__ == "__main__":
    test()