Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Compute bin_uiui via Flint for small n #383

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Compute bin_uiui via Flint for small n #383

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
41adb48
Compute bin_uiui via Flint for small n
albinahlback Oct 18, 2021
e811d3b
Use GMP's binomial in arb_bin_uiui
albinahlback Feb 26, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
109 changes: 87 additions & 22 deletions arb/bin.c
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
/*
Copyright (C) 2013 Fredrik Johansson
Copyright (C) 2021, 2022 Albin Ahlbäck

This file is part of Arb.

Expand All @@ -10,42 +11,106 @@
*/

#include "arb.h"
#include "flint/fmpz.h"

#define LOG2_BIN(n, k) \
(-1.3257480647361592 \
+ ((n) + .5) * log2(n) \
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think log2 might not be available everywhere. Better to use log.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

https://en.cppreference.com/w/c/numeric/math/log2

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"C99"

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is your opinion about supporting/assuming newer standards?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Might be needed at some point, but not worth risking any breakage right here.

- ((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);
}
else
_arb_bin_ui(x, n, k, prec);
}

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);
k = FLINT_MIN(k, n - k);

#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};
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This level of code shouldn't be be written in a way that relies on the internal representation of structs.

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
{
arf_struct atmp;
mag_struct mtmp = {0, 0};
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This level of code shouldn't be be written in a way that relies on the internal representation of structs.

arb_struct tmp = {atmp, mtmp};
arf_init_set_ui(&atmp, n);
arb_bin_ui(x, &tmp, k, prec);
}
}

#undef LOG2_BIN
8 changes: 6 additions & 2 deletions doc/source/arb.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down