From 41adb48f5abf2f874cdfa138960693ef168dbfa0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Albin=20Ahlb=C3=A4ck?= Date: Mon, 18 Oct 2021 17:39:41 +0200 Subject: [PATCH 1/2] Compute bin_uiui via Flint for small n --- arb/bin.c | 27 +++++++++++++++++++++++---- doc/source/arb.rst | 8 ++++++-- 2 files changed, 29 insertions(+), 6 deletions(-) diff --git a/arb/bin.c b/arb/bin.c index dc5d0cdb3..000b21e85 100644 --- a/arb/bin.c +++ b/arb/bin.c @@ -1,5 +1,6 @@ /* Copyright (C) 2013 Fredrik Johansson + Copyright (C) 2021 Albin Ahlbäck This file is part of Arb. @@ -10,6 +11,7 @@ */ #include "arb.h" +#include "flint/fmpz.h" void arb_bin_ui(arb_t x, const arb_t n, ulong k, slong prec) @@ -39,13 +41,30 @@ arb_bin_ui(arb_t x, const arb_t n, ulong k, slong prec) } } +static +void +_arb_bin_uiui_small(arb_t x, ulong n, ulong k, slong prec) +{ + fmpz_t b; + fmpz_init(b); + fmpz_bin_uiui(b, n, k); + arb_set_round_fmpz(x, b, prec); + fmpz_clear(b); +} + void arb_bin_uiui(arb_t x, ulong n, ulong k, slong prec) { arb_t t; - arb_init(t); - arb_set_ui(t, n); - arb_bin_ui(x, t, k, prec); - arb_clear(t); + + if (n < 3000) + _arb_bin_uiui_small(x, n, k, prec); + else + { + arb_init(t); + arb_set_ui(t, n); + arb_bin_ui(x, t, k, prec); + arb_clear(t); + } } diff --git a/doc/source/arb.rst b/doc/source/arb.rst index 039882bf9..9a4e09a34 100644 --- a/doc/source/arb.rst +++ b/doc/source/arb.rst @@ -1348,11 +1348,15 @@ Gamma function and factorials .. function:: void arb_bin_ui(arb_t z, const arb_t n, ulong k, slong prec) -.. function:: void arb_bin_uiui(arb_t z, ulong n, ulong k, slong prec) - Computes the binomial coefficient `z = {n \choose k}`, via the rising factorial as `{n \choose k} = (n-k+1)_k / k!`. +.. function:: void arb_bin_uiui(arb_t z, ulong n, ulong k, slong prec) + + Computes the binomial coefficient `z = {n \choose k}`. If `n` is small, + it is computed via Flint's implementation `fmpz_bin_uiui`. If `n` is big, + it is computed via rising factorial as `{n \choose k} = (n-k+1)_k / k!`. + .. function:: void arb_gamma(arb_t z, const arb_t x, slong prec) void arb_gamma_fmpq(arb_t z, const fmpq_t x, slong prec) void arb_gamma_fmpz(arb_t z, const fmpz_t x, slong prec) From e811d3b6113abdc8b54b3784ccaf9f479185a967 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Albin=20Ahlb=C3=A4ck?= Date: Sat, 26 Feb 2022 14:46:41 +0100 Subject: [PATCH 2/2] Use GMP's binomial in arb_bin_uiui --- arb/bin.c | 118 +++++++++++++++++++++++++++++++++++++----------------- 1 file changed, 82 insertions(+), 36 deletions(-) diff --git a/arb/bin.c b/arb/bin.c index 000b21e85..d75284683 100644 --- a/arb/bin.c +++ b/arb/bin.c @@ -1,6 +1,6 @@ /* Copyright (C) 2013 Fredrik Johansson - Copyright (C) 2021 Albin Ahlbäck + Copyright (C) 2021, 2022 Albin Ahlbäck This file is part of Arb. @@ -13,58 +13,104 @@ #include "arb.h" #include "flint/fmpz.h" +#define LOG2_BIN(n, k) \ + (-1.3257480647361592 \ + + ((n) + .5) * log2(n) \ + - ((k) + .5) * log2(k) \ + - ((n) - (k) + .5) * log2((n) - (k))) \ + +static void +_arb_bin_ui(arb_t x, const arb_t n, ulong k, slong prec) +{ + arb_t t, u; + + arb_init(t); + arb_init(u); + + arb_sub_ui(t, n, k - 1, prec); + arb_rising_ui(t, t, k, prec); + arb_fac_ui(u, k, prec); + arb_div(x, t, u, prec); + + arb_clear(t); + arb_clear(u); +} + void arb_bin_ui(arb_t x, const arb_t n, ulong k, slong prec) { if (k == 0) - { arb_one(x); - } else if (k == 1) - { arb_set_round(x, n, prec); - } - else + else if (arb_is_int(n) + && arb_is_nonnegative(n) + && ARF_EXP(arb_midref(n)) <= FLINT_BITS) /* fits in ulong */ { - arb_t t, u; - - arb_init(t); - arb_init(u); - - arb_sub_ui(t, n, k - 1, prec); - arb_rising_ui(t, t, k, prec); - arb_fac_ui(u, k, prec); - arb_div(x, t, u, prec); - - arb_clear(t); - arb_clear(u); + fmpz_t tmp; + fmpz_init(tmp); + arb_get_unique_fmpz(tmp, n); + arb_bin_uiui(x, fmpz_get_ui(tmp), k, prec); } -} - -static -void -_arb_bin_uiui_small(arb_t x, ulong n, ulong k, slong prec) -{ - fmpz_t b; - fmpz_init(b); - fmpz_bin_uiui(b, n, k); - arb_set_round_fmpz(x, b, prec); - fmpz_clear(b); + else + _arb_bin_ui(x, n, k, prec); } void arb_bin_uiui(arb_t x, ulong n, ulong k, slong prec) { - arb_t t; + k = FLINT_MIN(k, n - k); - if (n < 3000) - _arb_bin_uiui_small(x, n, k, prec); +#if LONG_MAX == WORD_MAX + if (n <= (UWORD(1) << 12) + || k <= (UWORD(1) << 7) +#if FLINT_BITS == 64 + || (n <= (UWORD(1) << 32) && k <= (UWORD(1) << 8)) +#else + || (n <= UWORD_MAX && k <= (UWORD(1) << 8)) +#endif + || (n <= (UWORD(1) << 22) && k <= (UWORD(1) << 9)) + || (n <= (UWORD(1) << 16) && k <= (UWORD(1) << 10)) + || ((double) prec) < LOG2_BIN(n, k)) + { + mpz_t mres; + mpz_init(mres); + flint_mpz_bin_uiui(mres, n, k); + arb_set_round_mpz(x, mres, prec); + mpz_clear(mres); + } +#else + if (n <= (UWORD(1) << 12) + || k <= (UWORD(1) << 7) + || (n < UWORD_MAX && k <= (UWORD(1) << 8)) + || (n <= (UWORD(1) << 22) && k <= (UWORD(1) << 9)) + || (n <= (UWORD(1) << 16) && k <= (UWORD(1) << 10))) + { + mpz_t mres; + mpz_init(mres); + flint_mpz_bin_uiui(mres, n, k); + arb_set_round_mpz(x, mres, prec); + mpz_clear(mres); + } + else if (k < (UWORD(1) << 32) && ((double) prec) < LOG2_BIN(n, k)) + { + mpz_t mres; + __mpz_struct mn = {1, 1, NULL}; + mpz_init(mres); + mn._mp_d = &n; + mpz_bin_ui(mres, &mn, k); + arb_set_round_mpz(x, mres, prec); + mpz_clear(mres); + } +#endif else { - arb_init(t); - arb_set_ui(t, n); - arb_bin_ui(x, t, k, prec); - arb_clear(t); + arf_struct atmp; + mag_struct mtmp = {0, 0}; + arb_struct tmp = {atmp, mtmp}; + arf_init_set_ui(&atmp, n); + arb_bin_ui(x, &tmp, k, prec); } } +#undef LOG2_BIN