Binomial coefficient of two unsigned long long values avoiding overflow if possible

binomial.c

#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <limits.h>

/**
 * \brief Returns Greatest Common Divisor of two unsigned long long values.
 *
 * \param[in] a First value.
 * \param[in] b Second value.
 * 
 * \returns GCD of \p a and \p b if both are non-zero.
 * \returns \p a if \p b is zero.
 * \returns \p b if \p a is zero.
 */
unsigned long long gcd_ull(unsigned long long a, unsigned long long b)
{
    while (b) {
        unsigned long long t = b;

        b = a % b;
        a = t;
    }
    return a;
}

/**
 * \brief Calculates binomial coefficient of two unsigned long long values.
 *
 * Calculates a binomial coefficient as long as the arguments are valid
 * and the result is no more than \c ULLONG_MAX.
 *
 * \param[in] n
 * \param[in] k
 * \param[out] res Pointer to storage for result.
 *
 * \returns 0 on success.
 * \returns -1 on error and sets \c errno.
 * 
 * \par Errors
 *
 * \c errno  | Cause
 * --------  | -----
 * \c EINVAL | Invalid arguments (\p k > \p n).
 * \c ERANGE | Result is out of range.
 */
int calc_binomial_ull(unsigned long long n, unsigned long long k,
                      unsigned long long *res)
{
    unsigned long long v;
    unsigned long long i;

    *res = 0;
    if (k > n) {
        errno = EINVAL;
        return -1;
    }
    if (k > n / 2) {
        k = n - k;
    }
    v = 1;
    for (i = 1; i <= k; i++) {
        unsigned long long num = n - k + i;
        unsigned long long den = i;

        if (ULLONG_MAX / num < v) {
            /*
             * Try and avoid overflow by reducing numerator and
             * denominator to smallest terms and doing the division
             * before the multiplication.
             */
            unsigned long long g = gcd_ull(num, den);

            num /= g;
            den /= g;
            v /= den;
            if (ULLONG_MAX / num < v) {
                /* It will still overflow. */
                errno = ERANGE;
                return -1;
            }
            v *= num;
        } else {
            v = v * num / den;
        }
    }
    *res = v;
    return 0;
}

static const char *progname = "binomial";

static void usage(void)
{
    fprintf(stderr,
        "usage: %s N K\n"
        "Calculate N choose K\n", progname);
}

static int argull(const char *arg, int base, unsigned long long *val)
{
    char *endptr;

    errno = 0;
    *val = strtoull(arg, &endptr, base);
    if (endptr == arg || *endptr != '\0') {
        errno = EINVAL;
    }
    if (errno) {
        if (errno == ERANGE) {
            return -1;
        }
        return -2;
    }
    return 0;
}

int main(int argc, char *argv[])
{
    unsigned long long n;
    unsigned long long k;
    unsigned long long nck;
    int rc = 2;

    if (argc) {
        progname = argv[0];
    }
    if (argc != 3) {
        usage();
        return rc;
    }
    rc = -argull(argv[1], 10, &n);
    if (rc) {
        perror(argv[1]);
        if (rc == 2) {
            usage();
        }
        return rc;
    }
    rc = -argull(argv[2], 10, &k);
    if (rc) {
        perror(argv[2]);
        if (rc == 2) {
            usage();
        }
        return rc;
    }
    rc = -calc_binomial_ull(n, k, &nck);
    if (rc) {
        perror(NULL);
        return rc;
    }
    printf("%llu\n", nck);
    return rc;
}