incomplete beta function

This commit is contained in:
Fredrik Johansson 2016-05-19 15:43:15 +02:00
parent f0dbadc14b
commit 1835534cae
6 changed files with 524 additions and 4 deletions

View file

@ -148,6 +148,10 @@ void acb_hypgeom_gamma_lower(acb_t res, const acb_t s, const acb_t z, int modifi
void _acb_hypgeom_gamma_lower_series(acb_ptr g, const acb_t s, acb_srcptr h, slong hlen, int regularized, slong n, slong prec);
void acb_hypgeom_gamma_lower_series(acb_poly_t g, const acb_t s, const acb_poly_t h, int regularized, slong n, slong prec);
void acb_hypgeom_beta_lower(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec);
void _acb_hypgeom_beta_lower_series(acb_ptr res, const acb_t a, const acb_t b, acb_srcptr z, slong zlen, int regularized, slong len, slong prec);
void acb_hypgeom_beta_lower_series(acb_poly_t res, const acb_t a, const acb_t b, const acb_poly_t z, int regularized, slong len, slong prec);
void acb_hypgeom_expint(acb_t res, const acb_t s, const acb_t z, slong prec);
void acb_hypgeom_erf_propagated_error(mag_t re, mag_t im, const acb_t z);

86
acb_hypgeom/beta_lower.c Normal file
View file

@ -0,0 +1,86 @@
/*
Copyright (C) 2016 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 "acb_hypgeom.h"
/* todo: move this to the acb module? */
static void
acb_beta(acb_t res, const acb_t a, const acb_t b, slong prec)
{
acb_t t, u;
acb_init(t);
acb_init(u);
acb_gamma(t, a, prec);
acb_gamma(u, b, prec);
acb_add(res, a, b, prec);
acb_rgamma(res, res, prec);
acb_mul(res, res, t, prec);
acb_mul(res, res, u, prec);
acb_clear(t);
acb_clear(u);
}
void acb_hypgeom_beta_lower(acb_t res,
const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec)
{
acb_t t, u;
if (acb_is_zero(z) && arb_is_positive(acb_realref(a)))
{
acb_zero(res);
return;
}
if (acb_is_one(z) && arb_is_positive(acb_realref(b)))
{
if (regularized)
acb_one(res);
else
acb_beta(res, a, b, prec);
return;
}
acb_init(t);
acb_init(u);
acb_sub_ui(t, b, 1, prec);
acb_neg(t, t);
acb_add_ui(u, a, 1, prec);
if (regularized)
{
acb_hypgeom_2f1(t, a, t, u, z, 1, prec);
acb_add(u, a, b, prec);
acb_gamma(u, u, prec);
acb_mul(t, t, u, prec);
acb_rgamma(u, b, prec);
acb_mul(t, t, u, prec);
}
else
{
acb_hypgeom_2f1(t, a, t, u, z, 0, prec);
acb_div(t, t, a, prec);
}
acb_pow(u, z, a, prec);
acb_mul(t, t, u, prec);
acb_set(res, t);
acb_clear(t);
acb_clear(u);
}

View file

@ -0,0 +1,114 @@
/*
Copyright (C) 2016 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 "acb_hypgeom.h"
void
_acb_hypgeom_beta_lower_series(acb_ptr res,
const acb_t a, const acb_t b, acb_srcptr z, slong zlen, int regularized,
slong len, slong prec)
{
acb_ptr t, u, v;
acb_t c, d, e;
zlen = FLINT_MIN(zlen, len);
if (zlen == 1)
{
acb_hypgeom_beta_lower(res, a, b, z, regularized, prec);
_acb_vec_zero(res + 1, len - 1);
return;
}
t = _acb_vec_init(len);
u = _acb_vec_init(len);
v = _acb_vec_init(zlen - 1);
acb_init(c);
acb_init(d);
acb_init(e);
acb_hypgeom_beta_lower(d, a, b, z, regularized, prec);
if (regularized)
{
/* todo: except in special cases, we already compute a bunch of
gamma functions in beta_lower, so we could avoid recomputing them */
acb_add(e, a, b, prec);
acb_gamma(e, e, prec);
acb_rgamma(c, a, prec);
acb_mul(e, e, c, prec);
acb_rgamma(c, b, prec);
acb_mul(e, e, c, prec);
}
/* u = (1-z)^(b-1) */
_acb_vec_neg(t, z, zlen);
acb_add_ui(t, t, 1, prec);
acb_sub_ui(c, b, 1, prec);
_acb_poly_pow_acb_series(u, t, FLINT_MIN(zlen, len - 1), c, len - 1, prec);
/* t = z^(a-1) */
acb_sub_ui(c, a, 1, prec);
_acb_poly_pow_acb_series(t, z, FLINT_MIN(zlen, len - 1), c, len - 1, prec);
/* v = z' */
_acb_poly_derivative(v, z, zlen, prec);
_acb_poly_mullow(res, t, len - 1, u, len - 1, len - 1, prec);
_acb_poly_mullow(t, res, len - 1, v, zlen - 1, len - 1, prec);
_acb_poly_integral(res, t, len, prec);
if (regularized)
{
_acb_vec_scalar_mul(res, res, len, e, prec);
}
acb_set(res, d);
_acb_vec_clear(t, len);
_acb_vec_clear(u, len);
_acb_vec_clear(v, zlen - 1);
acb_clear(c);
acb_clear(d);
acb_clear(e);
}
void acb_hypgeom_beta_lower_series(acb_poly_t res,
const acb_t a, const acb_t b, const acb_poly_t z, int regularized,
slong len, slong prec)
{
if (len == 0)
{
acb_poly_zero(res);
return;
}
acb_poly_fit_length(res, len);
if (z->length == 0)
{
acb_t t;
acb_init(t);
_acb_hypgeom_beta_lower_series(res->coeffs, a, b, t, 1,
regularized, len, prec);
acb_clear(t);
}
else
{
_acb_hypgeom_beta_lower_series(res->coeffs, a, b, z->coeffs,
z->length, regularized, len, prec);
}
_acb_poly_set_length(res, len);
_acb_poly_normalise(res);
}

View file

@ -0,0 +1,130 @@
/*
Copyright (C) 2016 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 "acb_hypgeom.h"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("beta_lower....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
{
acb_t a, a1, b, b1, z, w, wa1, wb1, t, u;
slong prec;
int regularized;
acb_init(a);
acb_init(a1);
acb_init(b);
acb_init(b1);
acb_init(z);
acb_init(w);
acb_init(wa1);
acb_init(wb1);
acb_init(t);
acb_init(u);
regularized = n_randint(state, 2);
prec = 2 + n_randint(state, 100);
acb_randtest_param(a, state, 1 + n_randint(state, 100), 1 + n_randint(state, 10));
acb_randtest_param(b, state, 1 + n_randint(state, 100), 1 + n_randint(state, 10));
acb_randtest_param(z, state, 1 + n_randint(state, 100), 1 + n_randint(state, 10));
acb_add_ui(a1, a, 1, prec);
acb_add_ui(b1, b, 1, prec);
acb_hypgeom_beta_lower(w, a, b, z, regularized, prec);
acb_hypgeom_beta_lower(wa1, a1, b, z, regularized, prec);
acb_hypgeom_beta_lower(wb1, a, b1, z, regularized, prec);
if (regularized)
{
acb_mul(t, wa1, a, prec);
acb_addmul(t, wb1, b, prec);
acb_add(u, a, b, prec);
acb_div(t, t, u, prec);
}
else
{
acb_add(t, wa1, wb1, prec);
}
if (!acb_overlaps(w, t))
{
flint_printf("FAIL: contiguous relation\n\n");
flint_printf("regularized = %d\n\n", regularized);
flint_printf("a = "); acb_printd(a, 30); flint_printf("\n\n");
flint_printf("b = "); acb_printd(b, 30); flint_printf("\n\n");
flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n");
flint_printf("w = "); acb_printd(w, 30); flint_printf("\n\n");
flint_printf("wa1 = "); acb_printd(wa1, 30); flint_printf("\n\n");
flint_printf("wb1 = "); acb_printd(wb1, 30); flint_printf("\n\n");
flint_printf("t = "); acb_printd(t, 30); flint_printf("\n\n");
abort();
}
/* test I(a,b;z) = 1-I(b,a,1-z) */
if (!regularized)
{
acb_add(t, a, b, prec);
acb_gamma(t, t, prec);
acb_mul(w, w, t, prec);
acb_rgamma(t, a, prec);
acb_mul(w, w, t, prec);
acb_rgamma(t, b, prec);
acb_mul(w, w, t, prec);
}
acb_sub_ui(t, z, 1, prec);
acb_neg(t, t);
acb_hypgeom_beta_lower(t, b, a, t, 1, prec);
acb_sub_ui(t, t, 1, prec);
acb_neg(t, t);
if (!acb_overlaps(w, t))
{
flint_printf("FAIL: symmetry\n\n");
flint_printf("regularized = %d\n\n", regularized);
flint_printf("a = "); acb_printd(a, 30); flint_printf("\n\n");
flint_printf("b = "); acb_printd(b, 30); flint_printf("\n\n");
flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n");
flint_printf("w = "); acb_printd(w, 30); flint_printf("\n\n");
flint_printf("t = "); acb_printd(t, 30); flint_printf("\n\n");
abort();
}
acb_clear(a);
acb_clear(a1);
acb_clear(b);
acb_clear(b1);
acb_clear(z);
acb_clear(w);
acb_clear(wa1);
acb_clear(wb1);
acb_clear(t);
acb_clear(u);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

View file

@ -0,0 +1,155 @@
/*
Copyright (C) 2016 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 "acb_hypgeom.h"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("beta_lower_series....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
{
acb_t a, a1, b, b1, t, u;
acb_poly_t z, w, wa1, wb1, pt, pu;
slong prec, len, lena1, lenb1;
int regularized;
acb_init(a);
acb_init(a1);
acb_init(b);
acb_init(b1);
acb_init(t);
acb_init(u);
acb_poly_init(z);
acb_poly_init(w);
acb_poly_init(wa1);
acb_poly_init(wb1);
acb_poly_init(pt);
acb_poly_init(pu);
regularized = n_randint(state, 2);
prec = 2 + n_randint(state, 100);
acb_randtest_param(a, state, 1 + n_randint(state, 100), 1 + n_randint(state, 10));
acb_randtest_param(b, state, 1 + n_randint(state, 100), 1 + n_randint(state, 10));
acb_poly_randtest(z, state, 10, 1 + n_randint(state, 100), 1 + n_randint(state, 10));
acb_poly_randtest(w, state, 10, 1 + n_randint(state, 100), 1 + n_randint(state, 10));
acb_poly_randtest(wa1, state, 10, 1 + n_randint(state, 100), 1 + n_randint(state, 10));
acb_poly_randtest(wb1, state, 10, 1 + n_randint(state, 100), 1 + n_randint(state, 10));
len = n_randint(state, 10);
lena1 = n_randint(state, 10);
lenb1 = n_randint(state, 10);
acb_add_ui(a1, a, 1, prec);
acb_add_ui(b1, b, 1, prec);
acb_hypgeom_beta_lower_series(w, a, b, z, regularized, len, prec);
acb_hypgeom_beta_lower_series(wa1, a1, b, z, regularized, lena1, prec);
acb_hypgeom_beta_lower_series(wb1, a, b1, z, regularized, lenb1, prec);
if (regularized)
{
acb_poly_scalar_mul(pt, wa1, a, prec);
acb_poly_scalar_mul(pu, wb1, b, prec);
acb_poly_add(pt, pt, pu, prec);
acb_add(u, a, b, prec);
acb_poly_scalar_div(pt, pt, u, prec);
}
else
{
acb_poly_add(pt, wa1, wb1, prec);
}
acb_poly_set(pu, w);
acb_poly_truncate(pu, FLINT_MIN(FLINT_MIN(len, lena1), lenb1));
acb_poly_truncate(pt, FLINT_MIN(FLINT_MIN(len, lena1), lenb1));
if (!acb_poly_overlaps(pu, pt))
{
flint_printf("FAIL: contiguous relation\n\n");
flint_printf("regularized = %d\n\n", regularized);
flint_printf("len = %wd, lena1 = %wd, lenb1 = %wd\n\n", len, lena1, lenb1);
flint_printf("a = "); acb_printd(a, 30); flint_printf("\n\n");
flint_printf("b = "); acb_printd(b, 30); flint_printf("\n\n");
flint_printf("z = "); acb_poly_printd(z, 30); flint_printf("\n\n");
flint_printf("w = "); acb_poly_printd(w, 30); flint_printf("\n\n");
flint_printf("wa1 = "); acb_poly_printd(wa1, 30); flint_printf("\n\n");
flint_printf("wb1 = "); acb_poly_printd(wb1, 30); flint_printf("\n\n");
flint_printf("pt = "); acb_poly_printd(pt, 30); flint_printf("\n\n");
abort();
}
/* test I(a,b;z) = 1-I(b,a,1-z) */
if (!regularized)
{
acb_add(t, a, b, prec);
acb_gamma(t, t, prec);
acb_poly_scalar_mul(w, w, t, prec);
acb_rgamma(t, a, prec);
acb_poly_scalar_mul(w, w, t, prec);
acb_rgamma(t, b, prec);
acb_poly_scalar_mul(w, w, t, prec);
}
acb_poly_add_si(pt, z, -1, prec);
acb_poly_neg(pt, pt);
acb_hypgeom_beta_lower_series(pt, b, a, pt, 1, lena1, prec);
acb_poly_add_si(pt, pt, -1, prec);
acb_poly_neg(pt, pt);
acb_poly_set(pu, w);
acb_poly_truncate(pu, FLINT_MIN(len, lena1));
acb_poly_truncate(pt, FLINT_MIN(len, lena1));
if (!acb_poly_overlaps(pu, pt))
{
flint_printf("FAIL: symmetry\n\n");
flint_printf("regularized = %d\n\n", regularized);
flint_printf("len = %wd, lena1 = %wd\n\n", len, lena1);
flint_printf("a = "); acb_printd(a, 30); flint_printf("\n\n");
flint_printf("b = "); acb_printd(b, 30); flint_printf("\n\n");
flint_printf("z = "); acb_poly_printd(z, 30); flint_printf("\n\n");
flint_printf("w = "); acb_poly_printd(w, 30); flint_printf("\n\n");
flint_printf("pt = "); acb_poly_printd(pt, 30); flint_printf("\n\n");
abort();
}
acb_clear(a);
acb_clear(a1);
acb_clear(b);
acb_clear(b1);
acb_clear(t);
acb_clear(u);
acb_poly_clear(z);
acb_poly_clear(w);
acb_poly_clear(wa1);
acb_poly_clear(wb1);
acb_poly_clear(pt);
acb_poly_clear(pu);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

View file

@ -560,7 +560,7 @@ simultaneously. Any of the four function values can be omitted by passing
truncated to length *len*. As with the other Airy methods, any of the
outputs can be *NULL*.
Incomplete gamma functions
Incomplete gamma and beta functions
-------------------------------------------------------------------------------
.. function:: void acb_hypgeom_gamma_upper_asymp(acb_t res, const acb_t s, const acb_t z, int regularized, slong prec)
@ -610,7 +610,7 @@ Incomplete gamma functions
The *singular* version evaluates the finite sum directly and therefore
assumes that *s* is not too large.
.. function:: void _acb_hypgeom_gamma_upper_series(acb_ptr res, acb_t s, acb_srcptr z, slong zlen, int regularized, slong n, slong prec)
.. function:: void _acb_hypgeom_gamma_upper_series(acb_ptr res, const acb_t s, acb_srcptr z, slong zlen, int regularized, slong n, slong prec)
.. function:: void acb_hypgeom_gamma_upper_series(acb_poly_t res, const acb_t s, const acb_poly_t z, int regularized, slong n, slong prec)
@ -619,7 +619,6 @@ Incomplete gamma functions
The *regularized* argument has the same interpretation as in
:func:`acb_hypgeom_gamma_upper`.
.. function:: void acb_hypgeom_gamma_lower(acb_t res, const acb_t s, const acb_t z, int regularized, slong prec)
If *regularized* is 0, computes the lower incomplete gamma function
@ -631,7 +630,7 @@ Incomplete gamma functions
If *regularized* is 2, computes a further regularized lower incomplete
gamma function `\gamma^{*}(s,z) = z^{-s} P(s,z)`.
.. function:: void _acb_hypgeom_gamma_lower_series(acb_ptr res, acb_t s, acb_srcptr z, slong zlen, int regularized, slong n, slong prec)
.. function:: void _acb_hypgeom_gamma_lower_series(acb_ptr res, const acb_t s, acb_srcptr z, slong zlen, int regularized, slong n, slong prec)
.. function:: void acb_hypgeom_gamma_lower_series(acb_poly_t res, const acb_t s, const acb_poly_t z, int regularized, slong n, slong prec)
@ -640,6 +639,38 @@ Incomplete gamma functions
The *regularized* argument has the same interpretation as in
:func:`acb_hypgeom_gamma_lower`.
.. function:: void acb_hypgeom_beta_lower(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec)
Computes the (lower) incomplete beta function, defined by
`B(a,b;z) = \int_0^z t^{a-1} (1-t)^{b-1}`,
optionally the regularized incomplete beta function
`I(a,b;z) = B(a,b;z) / B(a,b;1)`.
In general, the integral must be interpreted using analytic continuation.
The precise definitions for all parameter values are
.. math ::
B(a,b;z) = \frac{z^a}{a} {}_2F_1(a, 1-b, a+1, z)
.. math ::
I(a,b;z) = \frac{\Gamma(a+b)}{\Gamma(b)} z^a {}_2{\widetilde F}_1(a, 1-b, a+1, z).
Note that both functions with this definition are undefined
for nonpositive integer *a*, and *I* is undefined for nonpositive integer
`a + b`.
.. function:: void _acb_hypgeom_beta_lower_series(acb_ptr res, const acb_t a, const acb_t b, acb_srcptr z, slong zlen, int regularized, slong n, slong prec)
.. function:: void acb_hypgeom_beta_lower_series(acb_poly_t res, const acb_t a, const acb_t b, const acb_poly_t z, int regularized, slong n, slong prec)
Sets *res* to the lower incomplete beta function `B(a,b;z)` (optionally
the regularized version `I(a,b;z)`) where *a* and *b* are constants
and *z* is a power series, truncating the result to length *n*.
The underscore method requires positive lengths and does not support
aliasing.
Exponential and trigonometric integrals
-------------------------------------------------------------------------------