add mag_bin_uiui

doc/source/mag.rst

@@ -384,6 +384,10 @@ Special functions
 
     Sets *z* to an upper bound for `1/n!`.
 
+.. function:: void mag_bin_uiui(mag_t res, ulong n, ulong k)
+
+    Sets *res* to an upper bound for the binomial coefficient `{n \choose k}`.
+
 .. function:: void mag_bernoulli_div_fac_ui(mag_t z, ulong n)
 
     Sets *z* to an upper bound for `|B_n| / n!` where `B_n` denotes

mag.h

@@ -607,6 +607,7 @@ void mag_pow_fmpz(mag_t z, const mag_t x, const fmpz_t e);
 
 void mag_fac_ui(mag_t z, ulong n);
 void mag_rfac_ui(mag_t z, ulong n);
+void mag_bin_uiui(mag_t res, ulong n, ulong k);
 
 /* TODO: test */
 void mag_bernoulli_div_fac_ui(mag_t z, ulong n);

mag/bin_uiui.c

@@ -0,0 +1,69 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "mag.h"
+
+void
+mag_bin_uiui(mag_t res, ulong n, ulong k)
+{
+    mag_t t;
+
+    if (k > n)
+    {
+        mag_zero(res);
+        return;
+    }
+
+    if (k > n / 2)
+        k = n - k;
+
+    if (k == 0)
+    {
+        mag_one(res);
+        return;
+    }
+
+    if (k == 1)
+    {
+        mag_set_ui(res, n);
+        return;
+    }
+
+    mag_init(t);
+
+    if (n < 256 && k < 256)
+    {
+        /* using accurate (lookup table) factorials */
+        mag_fac_ui(res, n);
+        mag_rfac_ui(t, k);
+        mag_mul(res, res, t);
+        mag_rfac_ui(t, n - k);
+        mag_mul(res, res, t);
+    }
+    else
+    {
+        /* binary entropy bound (n/(n-k))^(n-k) (n/k)^k = n^n / (k^k (n-k)^(n-k)) */
+        mag_set_ui(res, n);
+        mag_div_ui(res, res, n - k);
+        mag_pow_ui(res, res, n - k);
+        /* (n/(n-k))^(n-k) is also bounded by exp(k) */
+        mag_set_ui(t, k);
+        mag_exp(t, t);
+        mag_min(res, res, t);
+        mag_set_ui(t, n);
+        mag_div_ui(t, t, k);
+        mag_pow_ui(t, t, k);
+        mag_mul(res, res, t);
+    }
+
+    mag_clear(t);
+}
+

mag/test/t-bin_uiui.c

@@ -0,0 +1,67 @@
+/*
+    Copyright (C) 2014 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "mag.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("bin_ui....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 2000 * arb_test_multiplier(); iter++)
+    {
+        fmpr_t x, y;
+        fmpz_t f;
+        mag_t xb;
+        ulong n, k;
+
+        fmpr_init(x);
+        fmpr_init(y);
+        fmpz_init(f);
+        mag_init(xb);
+
+        mag_randtest_special(xb, state, 80);
+        n = n_randtest(state) % 10000;
+        k = n_randtest(state) % 10000;
+
+        mag_bin_uiui(xb, n, k);
+        fmpz_bin_uiui(f, n, k);
+        fmpr_set_fmpz(x, f);
+        mag_get_fmpr(y, xb);
+
+        MAG_CHECK_BITS(xb)
+
+        if (!(fmpr_cmpabs(y, x) >= 0))
+        {
+            flint_printf("FAIL\n\n");
+            flint_printf("n = %wu\n\n", n);
+            flint_printf("x = "); fmpr_print(x); flint_printf("\n\n");
+            flint_printf("y = "); fmpr_print(y); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        fmpr_clear(x);
+        fmpr_clear(y);
+        fmpz_clear(f);
+        mag_clear(xb);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+