hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-06 01:41:39 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

refine fmprb test functions; add gamma/rgamma/lgamma with test code

Browse source
This commit is contained in:
Fredrik Johansson 2012-12-22 21:29:11 +01:00
parent 7a310917d7
commit 147297748f
7 changed files with 367 additions and 23 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

27
doc/source/fmprb.rst
View file

@ -178,6 +178,20 @@ Random number generation
Generates a random ball. The midpoint and radius will both be finite.
.. function:: void fmprb_randtest_exact(fmprb_t x, flint_rand_t state, long prec, long mag_bits)
Generates a random number with zero radius.
.. function:: void fmprb_randtest_precise(fmprb_t x, flint_rand_t state, long prec, long mag_bits)
Generates a random number with radius at most `2^{-\mathrm{prec}}`
the magnitude of the midpoint.
.. function:: void fmprb_randtest_wide(fmprb_t x, flint_rand_t state, long prec, long mag_bits)
Generates a random number with midpoint and radius chosen independently,
possibly giving a very large interval.
.. function:: void fmprb_get_rand_fmpq(fmpq_t q, flint_rand_t state, const fmprb_t x, long bits)
Sets *q* to a random rational number from the interval represented by *x*.
@ -723,9 +737,14 @@ Gamma function
be done here so that a working precision is selected which always is
sufficient and also nearly optimal.
.. function:: void fmprb_gamma_log(fmprb_t y, const fmprb_t x, long prec)
.. function:: void fmprb_gamma(fmprb_t y, const fmprb_t x, long prec)
Sets `y = \log \Gamma(x)`, assuming that `x > 0`.
.. function:: void fmprb_rgamma(fmprb_t y, const fmprb_t x, long prec)
.. function:: void fmprb_lgamma(fmprb_t y, const fmprb_t x, long prec)
Sets, respectively, `y = \Gamma(x)`, `y = 1/\Gamma(x)`,
`y = \log \Gamma(x)`.
For large `x`, uses Stirling's expansion
@ -744,8 +763,8 @@ Gamma function
If `x` is too small for the asymptotic expansion to give sufficient
accuracy directly, we translate to `x + r`
using the formula `\log \Gamma(x) = \log \Gamma(x+r) -
\log(x (x+1) (x+2) \cdots (x+r-1))`.
using the formula `\Gamma(x) = \Gamma(x+r) /
(x (x+1) (x+2) \cdots (x+r-1))`.
To obtain a remainder smaller than `2^{-b}`, we must choose an `r` such
that `x + r > \beta b`, where `\beta > \log(2) / (2 \pi) \approx 0.11`.

8
fmprb.h
View file

@ -316,7 +316,10 @@ void fmprb_zeta_ui_vec_odd(fmprb_struct * x, ulong start, long num, long prec);
void fmprb_zeta_ui_vec(fmprb_struct * x, ulong start, long num, long prec);
void fmprb_gamma_fmpq_karatsuba(fmprb_struct * v, const fmpq_t a, long num, long prec);
void fmprb_gamma_log(fmprb_t y, const fmprb_t x, long prec);
void fmprb_lgamma(fmprb_t y, const fmprb_t x, long prec);
void fmprb_rgamma(fmprb_t y, const fmprb_t x, long prec);
void fmprb_gamma(fmprb_t y, const fmprb_t x, long prec);
static __inline__ void
fmprb_print(const fmprb_t x)
@ -493,6 +496,9 @@ void fmprb_add_error_2exp_si(fmprb_t x, long err);
void fmprb_add_error(fmprb_t x, const fmprb_t error);
void fmprb_randtest(fmprb_t x, flint_rand_t state, long prec, long mag_bits);
void fmprb_randtest_exact(fmprb_t x, flint_rand_t state, long prec, long mag_bits);
void fmprb_randtest_wide(fmprb_t x, flint_rand_t state, long prec, long mag_bits);
void fmprb_randtest_precise(fmprb_t x, flint_rand_t state, long prec, long mag_bits);
void fmprb_get_rand_fmpq(fmpq_t q, flint_rand_t state, const fmprb_t x, long bits);

55
fmprb/gamma.c
View file

@ -159,7 +159,7 @@ fmprb_gamma_log_stirling(fmprb_t s, const fmprb_t z, long nterms, long prec)
}
void
fmprb_gamma_log(fmprb_t y, const fmprb_t x, long prec)
fmprb_lgamma(fmprb_t y, const fmprb_t x, long prec)
{
long r, n, wp;
fmprb_t t, u;
@ -183,3 +183,56 @@ fmprb_gamma_log(fmprb_t y, const fmprb_t x, long prec)
fmprb_clear(t);
fmprb_clear(u);
}
void
fmprb_gamma(fmprb_t y, const fmprb_t x, long prec)
{
long r, n, wp;
fmprb_t t, u;
wp = prec + FLINT_BIT_COUNT(prec);
r = stirling_choose_r(x, wp);
n = stirling_choose_nterms(x, r, wp);
/* gamma(x) = gamma(x+r) / rf(x,r) */
fmprb_init(t);
fmprb_init(u);
fmprb_add_ui(t, x, r, wp);
fmprb_gamma_log_stirling(u, t, n, wp);
fmprb_exp(u, u, prec);
fmprb_rfac_ui_bsplit(t, x, r, wp);
fmprb_div(y, u, t, prec);
fmprb_clear(t);
fmprb_clear(u);
}
void
fmprb_rgamma(fmprb_t y, const fmprb_t x, long prec)
{
long r, n, wp;
fmprb_t t, u;
wp = prec + FLINT_BIT_COUNT(prec);
r = stirling_choose_r(x, wp);
n = stirling_choose_nterms(x, r, wp);
/* 1/gamma(x) = rf(x,r) / gamma(x+r) */
fmprb_init(t);
fmprb_init(u);
fmprb_add_ui(t, x, r, wp);
fmprb_gamma_log_stirling(u, t, n, wp);
fmprb_neg(u, u);
fmprb_exp(u, u, prec);
fmprb_rfac_ui_bsplit(t, x, r, wp);
fmprb_mul(y, u, t, prec);
fmprb_clear(t);
fmprb_clear(u);
}

57
fmprb/randtest.c
View file

@ -26,26 +26,49 @@
#include "fmprb.h"
void
fmprb_randtest(fmprb_t x, flint_rand_t state, long prec, long mag_bits)
fmprb_randtest_exact(fmprb_t x, flint_rand_t state, long prec, long mag_bits)
{
fmpr_randtest(fmprb_midref(x), state, prec, mag_bits);
fmpr_zero(fmprb_radref(x));
}
void
fmprb_randtest_wide(fmprb_t x, flint_rand_t state, long prec, long mag_bits)
{
fmpr_randtest(fmprb_midref(x), state, prec, mag_bits);
fmpr_randtest(fmprb_radref(x), state, FMPRB_RAD_PREC, mag_bits);
fmpr_abs(fmprb_radref(x), fmprb_radref(x));
}
void
fmprb_randtest_precise(fmprb_t x, flint_rand_t state, long prec, long mag_bits)
{
fmpr_randtest(fmprb_midref(x), state, prec, mag_bits);
switch (n_randint(state, 8))
if (fmpr_is_zero(fmprb_midref(x)) || (n_randint(state, 8) == 0))
{
/* exact */
case 0:
fmpr_zero(fmprb_radref(x));
break;
/* arbitrary radius */
case 1:
fmpr_randtest(fmprb_radref(x), state, FMPRB_RAD_PREC, mag_bits);
fmpr_abs(fmprb_radref(x), fmprb_radref(x));
break;
/* "typical" radius */
default:
fmpr_randtest_not_zero(fmprb_radref(x), state, FMPRB_RAD_PREC, 4);
fmpr_abs(fmprb_radref(x), fmprb_radref(x));
fmpz_add(fmpr_expref(fmprb_radref(x)), fmpr_expref(fmprb_radref(x)),
fmpr_expref(fmprb_midref(x)));
fmpr_zero(fmprb_radref(x));
}
else
{
fmpz_randtest_not_zero(fmpr_manref(fmprb_radref(x)), state, FMPRB_RAD_PREC);
fmpz_abs(fmpr_manref(fmprb_radref(x)), fmpr_manref(fmprb_radref(x)));
fmpz_sub_ui(fmpr_expref(fmprb_radref(x)), fmpr_expref(fmprb_midref(x)), prec + FMPRB_RAD_PREC);
}
}
void
fmprb_randtest(fmprb_t x, flint_rand_t state, long prec, long mag_bits)
{
switch (n_randint(state, 8))
{
case 0:
fmprb_randtest_exact(x, state, prec, mag_bits);
break;
case 1:
fmprb_randtest_wide(x, state, prec, mag_bits);
break;
default:
fmprb_randtest_precise(x, state, prec, mag_bits);
}
}

81
fmprb/test/t-gamma.c Normal file
View file

@ -0,0 +1,81 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2012 Fredrik Johansson
******************************************************************************/
#include "fmprb.h"
int main()
{
long iter;
flint_rand_t state;
printf("gamma....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 3000; iter++)
{
fmprb_t a, b, c;
long prec1, prec2;
prec1 = 2 + n_randint(state, 1000);
prec2 = prec1 + 30;
fmprb_init(a);
fmprb_init(b);
fmprb_init(c);
fmprb_randtest_precise(a, state, 1 + n_randint(state, 1000), 3);
fmprb_gamma(b, a, prec1);
fmprb_gamma(c, a, prec2);
if (!fmprb_overlaps(b, c))
{
printf("FAIL: overlap\n\n");
printf("a = "); fmprb_print(a); printf("\n\n");
printf("b = "); fmprb_print(b); printf("\n\n");
printf("c = "); fmprb_print(c); printf("\n\n");
abort();
}
fmprb_gamma(a, a, prec2);
if (!fmprb_equal(a, c))
{
printf("FAIL: aliasing\n\n");
abort();
}
fmprb_clear(a);
fmprb_clear(b);
fmprb_clear(c);
}
flint_randclear(state);
_fmpz_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

81
fmprb/test/t-lgamma.c Normal file
View file

@ -0,0 +1,81 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2012 Fredrik Johansson
******************************************************************************/
#include "fmprb.h"
int main()
{
long iter;
flint_rand_t state;
printf("lgamma....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 3000; iter++)
{
fmprb_t a, b, c;
long prec1, prec2;
prec1 = 2 + n_randint(state, 1000);
prec2 = prec1 + 30;
fmprb_init(a);
fmprb_init(b);
fmprb_init(c);
fmprb_randtest_precise(a, state, 1 + n_randint(state, 1000), 3);
fmprb_lgamma(b, a, prec1);
fmprb_lgamma(c, a, prec2);
if (!fmprb_overlaps(b, c))
{
printf("FAIL: overlap\n\n");
printf("a = "); fmprb_print(a); printf("\n\n");
printf("b = "); fmprb_print(b); printf("\n\n");
printf("c = "); fmprb_print(c); printf("\n\n");
abort();
}
fmprb_lgamma(a, a, prec2);
if (!fmprb_equal(a, c))
{
printf("FAIL: aliasing\n\n");
abort();
}
fmprb_clear(a);
fmprb_clear(b);
fmprb_clear(c);
}
flint_randclear(state);
_fmpz_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

81
fmprb/test/t-rgamma.c Normal file
View file

@ -0,0 +1,81 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2012 Fredrik Johansson
******************************************************************************/
#include "fmprb.h"
int main()
{
long iter;
flint_rand_t state;
printf("rgamma....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 3000; iter++)
{
fmprb_t a, b, c;
long prec1, prec2;
prec1 = 2 + n_randint(state, 1000);
prec2 = prec1 + 30;
fmprb_init(a);
fmprb_init(b);
fmprb_init(c);
fmprb_randtest_precise(a, state, 1 + n_randint(state, 1000), 3);
fmprb_rgamma(b, a, prec1);
fmprb_rgamma(c, a, prec2);
if (!fmprb_overlaps(b, c))
{
printf("FAIL: overlap\n\n");
printf("a = "); fmprb_print(a); printf("\n\n");
printf("b = "); fmprb_print(b); printf("\n\n");
printf("c = "); fmprb_print(c); printf("\n\n");
abort();
}
fmprb_rgamma(a, a, prec2);
if (!fmprb_equal(a, c))
{
printf("FAIL: aliasing\n\n");
abort();
}
fmprb_clear(a);
fmprb_clear(b);
fmprb_clear(c);
}
flint_randclear(state);
_fmpz_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}