diff --git a/doc/source/mag.rst b/doc/source/mag.rst index 5a8a093a..ea3100b8 100644 --- a/doc/source/mag.rst +++ b/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 diff --git a/mag.h b/mag.h index 4768386f..72fd6b89 100644 --- a/mag.h +++ b/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); diff --git a/mag/bin_uiui.c b/mag/bin_uiui.c new file mode 100644 index 00000000..bce683ea --- /dev/null +++ b/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 . +*/ + +#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); +} + diff --git a/mag/test/t-bin_uiui.c b/mag/test/t-bin_uiui.c new file mode 100644 index 00000000..49e985ba --- /dev/null +++ b/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 . +*/ + +#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; +} +