mirror of
https://github.com/vale981/arb
synced 2025-03-05 09:21:38 -05:00
complete implementation of hypergeometric U
This commit is contained in:
Fredrik Johansson
parent
42dcdb4f74
commit
f3e6d7772b
5 changed files with 434 additions and 14 deletions
|
|
@ -64,44 +64,38 @@ void acb_hypgeom_pfq_series_direct(acb_poly_t res,
|
|||
void acb_hypgeom_u_asymp(acb_t res, const acb_t a, const acb_t b,
|
||||
const acb_t z, long n, long prec);
|
||||
|
||||
void acb_hypgeom_u_1f1_series(acb_poly_t res,
|
||||
const acb_poly_t a, const acb_poly_t b, const acb_poly_t z,
|
||||
long len, long prec);
|
||||
|
||||
void acb_hypgeom_u_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, long prec);
|
||||
void acb_hypgeom_u(acb_t res, const acb_t a, const acb_t b, const acb_t z, long prec);
|
||||
|
||||
int acb_hypgeom_u_use_asymp(const acb_t z, long prec);
|
||||
|
||||
void acb_hypgeom_erf_1f1a(acb_t res, const acb_t z, long prec);
|
||||
|
||||
void acb_hypgeom_erf_1f1b(acb_t res, const acb_t z, long prec);
|
||||
|
||||
void acb_hypgeom_erf_asymp(acb_t res, const acb_t z, long prec, long prec2);
|
||||
|
||||
void acb_hypgeom_erf(acb_t res, const acb_t z, long prec);
|
||||
|
||||
void acb_hypgeom_bessel_j_0f1(acb_t res, const acb_t nu, const acb_t z, long prec);
|
||||
|
||||
void acb_hypgeom_bessel_j_asymp(acb_t res, const acb_t nu, const acb_t z, long prec);
|
||||
|
||||
void acb_hypgeom_bessel_j(acb_t res, const acb_t nu, const acb_t z, long prec);
|
||||
|
||||
void acb_hypgeom_bessel_k_0f1(acb_t res, const acb_t nu, const acb_t z, long prec);
|
||||
|
||||
void acb_hypgeom_bessel_k_0f1_series(acb_poly_t res, const acb_poly_t n, const acb_poly_t z, long len, long prec);
|
||||
|
||||
void acb_hypgeom_bessel_k_asymp(acb_t res, const acb_t nu, const acb_t z, long prec);
|
||||
|
||||
void acb_hypgeom_bessel_k(acb_t res, const acb_t nu, const acb_t z, long prec);
|
||||
|
||||
void acb_hypgeom_gamma_upper_asymp(acb_t res, const acb_t s, const acb_t z, int modified, long prec);
|
||||
|
||||
void acb_hypgeom_gamma_upper_1f1a(acb_t res, const acb_t s, const acb_t z, int modified, long prec);
|
||||
|
||||
void acb_hypgeom_gamma_upper_1f1b(acb_t res, const acb_t s, const acb_t z, int modified, long prec);
|
||||
|
||||
void acb_hypgeom_gamma_upper_singular(acb_t res, long s, const acb_t z, int modified, long prec);
|
||||
|
||||
void acb_hypgeom_gamma_upper(acb_t res, const acb_t s, const acb_t z, int modified, long prec);
|
||||
|
||||
void acb_hypgeom_expint(acb_t res, const acb_t s, const acb_t z, long prec);
|
||||
|
||||
void acb_hypgeom_erfc(acb_t res, const acb_t z, long prec);
|
||||
|
||||
void acb_hypgeom_erfi(acb_t res, const acb_t z, long prec);
|
||||
|
||||
void acb_hypgeom_ei_asymp(acb_t res, const acb_t z, long prec);
|
||||
|
|
|
|||
|
|
@ -67,7 +67,7 @@ acb_hypgeom_si_asymp(acb_t res, const acb_t z, long prec)
|
|||
|
||||
if (arb_is_zero(acb_realref(z)))
|
||||
{
|
||||
/* the function is real */
|
||||
/* the function is imaginary */
|
||||
arb_zero(acb_realref(t));
|
||||
}
|
||||
else if (arb_is_positive(acb_realref(z)))
|
||||
|
|
|
|||
200
acb_hypgeom/test/t-u.c
Normal file
200
acb_hypgeom/test/t-u.c
Normal file
|
|
@ -0,0 +1,200 @@
|
|||
/*=============================================================================
|
||||
|
||||
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) 2014 Fredrik Johansson
|
||||
|
||||
******************************************************************************/
|
||||
|
||||
#include "acb_hypgeom.h"
|
||||
|
||||
void
|
||||
acb_randtest_maybe_half_int(acb_t x, flint_rand_t state, long prec, long size)
|
||||
{
|
||||
if (n_randint(state, 8) == 0)
|
||||
{
|
||||
fmpz_t t;
|
||||
fmpz_init(t);
|
||||
fmpz_randtest(t, state, 1 + n_randint(state, prec));
|
||||
arb_set_fmpz(acb_realref(x), t);
|
||||
arb_zero(acb_imagref(x));
|
||||
acb_mul_2exp_si(x, x, -1);
|
||||
fmpz_clear(t);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_randtest(x, state, prec, size);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
acb_hypgeom_u_asymp_proper(acb_t res, const acb_t a, const acb_t b, const acb_t z, long prec)
|
||||
{
|
||||
acb_t t;
|
||||
acb_init(t);
|
||||
acb_pow(t, z, a, prec);
|
||||
acb_hypgeom_u_asymp(res, a, b, z, -1, prec);
|
||||
acb_div(res, res, t, prec);
|
||||
acb_clear(t);
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
long iter;
|
||||
flint_rand_t state;
|
||||
|
||||
printf("u....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
||||
for (iter = 0; iter < 2000; iter++)
|
||||
{
|
||||
acb_t a0, a1, a2, b, z, w0, w1, w2, t, u;
|
||||
long prec0, prec1, prec2;
|
||||
|
||||
acb_init(a0);
|
||||
acb_init(a1);
|
||||
acb_init(a2);
|
||||
acb_init(b);
|
||||
acb_init(z);
|
||||
acb_init(w0);
|
||||
acb_init(w1);
|
||||
acb_init(w2);
|
||||
acb_init(t);
|
||||
acb_init(u);
|
||||
|
||||
prec0 = 2 + n_randint(state, 700);
|
||||
prec1 = 2 + n_randint(state, 700);
|
||||
prec2 = 2 + n_randint(state, 700);
|
||||
|
||||
acb_randtest_maybe_half_int(a0, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100));
|
||||
acb_randtest_maybe_half_int(b, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100));
|
||||
acb_randtest(z, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100));
|
||||
acb_randtest(w0, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100));
|
||||
acb_randtest(w1, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100));
|
||||
acb_randtest(w2, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100));
|
||||
|
||||
acb_add_ui(a1, a0, 1, prec0);
|
||||
acb_add_ui(a2, a0, 2, prec0);
|
||||
|
||||
switch (n_randint(state, 3))
|
||||
{
|
||||
case 0:
|
||||
acb_hypgeom_u_asymp_proper(w0, a0, b, z, prec0);
|
||||
break;
|
||||
case 1:
|
||||
acb_hypgeom_u_1f1(w0, a0, b, z, prec0);
|
||||
break;
|
||||
default:
|
||||
acb_hypgeom_u(w0, a0, b, z, prec0);
|
||||
}
|
||||
|
||||
switch (n_randint(state, 3))
|
||||
{
|
||||
case 0:
|
||||
acb_hypgeom_u_asymp_proper(w1, a0, b, z, prec1);
|
||||
break;
|
||||
case 1:
|
||||
acb_hypgeom_u_1f1(w1, a0, b, z, prec1);
|
||||
break;
|
||||
default:
|
||||
acb_hypgeom_u(w1, a0, b, z, prec1);
|
||||
}
|
||||
|
||||
if (!acb_overlaps(w0, w1))
|
||||
{
|
||||
printf("FAIL: consistency\n\n");
|
||||
printf("a = "); acb_printd(a0, 30); printf("\n\n");
|
||||
printf("b = "); acb_printd(b, 30); printf("\n\n");
|
||||
printf("z = "); acb_printd(z, 30); printf("\n\n");
|
||||
printf("w0 = "); acb_printd(w0, 30); printf("\n\n");
|
||||
printf("w1 = "); acb_printd(w1, 30); printf("\n\n");
|
||||
abort();
|
||||
}
|
||||
|
||||
switch (n_randint(state, 3))
|
||||
{
|
||||
case 0:
|
||||
acb_hypgeom_u_asymp_proper(w1, a1, b, z, prec1);
|
||||
break;
|
||||
case 1:
|
||||
acb_hypgeom_u_1f1(w1, a1, b, z, prec1);
|
||||
break;
|
||||
default:
|
||||
acb_hypgeom_u(w1, a1, b, z, prec1);
|
||||
}
|
||||
|
||||
switch (n_randint(state, 3))
|
||||
{
|
||||
case 0:
|
||||
acb_hypgeom_u_asymp_proper(w2, a2, b, z, prec2);
|
||||
break;
|
||||
case 1:
|
||||
acb_hypgeom_u_1f1(w2, a2, b, z, prec2);
|
||||
break;
|
||||
default:
|
||||
acb_hypgeom_u(w2, a2, b, z, prec2);
|
||||
}
|
||||
|
||||
acb_set(t, w0);
|
||||
|
||||
acb_mul_2exp_si(u, a0, 1);
|
||||
acb_sub(u, u, b, prec0);
|
||||
acb_add(u, u, z, prec0);
|
||||
acb_add_ui(u, u, 2, prec0);
|
||||
|
||||
acb_submul(t, w1, u, prec0);
|
||||
|
||||
acb_sub(u, a2, b, prec0);
|
||||
acb_mul(u, u, a1, prec0);
|
||||
acb_addmul(t, w2, u, prec0);
|
||||
|
||||
if (!acb_contains_zero(t))
|
||||
{
|
||||
printf("FAIL: contiguous relation\n\n");
|
||||
printf("a = "); acb_printd(a0, 30); printf("\n\n");
|
||||
printf("b = "); acb_printd(b, 30); printf("\n\n");
|
||||
printf("z = "); acb_printd(z, 30); printf("\n\n");
|
||||
printf("w0 = "); acb_printd(w0, 30); printf("\n\n");
|
||||
printf("w1 = "); acb_printd(w1, 30); printf("\n\n");
|
||||
printf("w2 = "); acb_printd(w2, 30); printf("\n\n");
|
||||
printf("t = "); acb_printd(t, 30); printf("\n\n");
|
||||
abort();
|
||||
}
|
||||
|
||||
acb_clear(a0);
|
||||
acb_clear(a1);
|
||||
acb_clear(a2);
|
||||
acb_clear(b);
|
||||
acb_clear(z);
|
||||
acb_clear(w0);
|
||||
acb_clear(w1);
|
||||
acb_clear(w2);
|
||||
acb_clear(t);
|
||||
acb_clear(u);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
||||
201
acb_hypgeom/u.c
Normal file
201
acb_hypgeom/u.c
Normal file
|
|
@ -0,0 +1,201 @@
|
|||
/*=============================================================================
|
||||
|
||||
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) 2015 Fredrik Johansson
|
||||
|
||||
******************************************************************************/
|
||||
|
||||
#include "acb_hypgeom.h"
|
||||
|
||||
void
|
||||
acb_hypgeom_u_1f1_series(acb_poly_t res,
|
||||
const acb_poly_t a, const acb_poly_t b, const acb_poly_t z,
|
||||
long len, long prec)
|
||||
{
|
||||
acb_poly_t s, u, A, B;
|
||||
acb_poly_struct aa[3];
|
||||
arb_t c;
|
||||
long wlen;
|
||||
int singular;
|
||||
|
||||
acb_poly_init(s);
|
||||
acb_poly_init(u);
|
||||
acb_poly_init(A);
|
||||
acb_poly_init(B);
|
||||
acb_poly_init(aa + 0);
|
||||
acb_poly_init(aa + 1);
|
||||
acb_poly_init(aa + 2);
|
||||
arb_init(c);
|
||||
|
||||
singular = (b->length == 0) || acb_is_int(b->coeffs);
|
||||
wlen = len + (singular != 0);
|
||||
|
||||
/* A = rgamma(a-b+1) * 1F~1(a,b,z) */
|
||||
acb_poly_sub(u, a, b, prec);
|
||||
acb_poly_add_si(u, u, 1, prec);
|
||||
acb_poly_rgamma_series(A, u, wlen, prec);
|
||||
|
||||
/* todo: handle a = 1 efficiently */
|
||||
acb_poly_set(aa, a);
|
||||
acb_poly_set(aa + 1, b);
|
||||
acb_poly_one(aa + 2);
|
||||
acb_hypgeom_pfq_series_direct(s, aa, 1, aa + 1, 2, z, -1, wlen, prec, 1);
|
||||
acb_poly_mullow(A, A, s, wlen, prec);
|
||||
|
||||
/* B = rgamma(a) * 1F~1(a-b+1,2-b,z) * z^(1-b) */
|
||||
acb_poly_set(aa, u);
|
||||
acb_poly_add_si(aa + 1, b, -2, prec);
|
||||
acb_poly_neg(aa + 1, aa + 1);
|
||||
acb_hypgeom_pfq_series_direct(s, aa, 1, aa + 1, 2, z, -1, wlen, prec, 1);
|
||||
acb_poly_rgamma_series(B, a, wlen, prec);
|
||||
acb_poly_mullow(B, B, s, wlen, prec);
|
||||
|
||||
acb_poly_add_si(u, b, -1, prec);
|
||||
acb_poly_neg(u, u);
|
||||
acb_poly_pow_series(s, z, u, wlen, prec);
|
||||
acb_poly_mullow(B, B, s, wlen, prec);
|
||||
|
||||
acb_poly_sub(A, A, B, prec);
|
||||
|
||||
/* multiply by pi csc(pi b) */
|
||||
acb_poly_sin_pi_series(B, b, wlen, prec);
|
||||
|
||||
if (singular)
|
||||
{
|
||||
acb_poly_shift_right(A, A, 1);
|
||||
acb_poly_shift_right(B, B, 1);
|
||||
}
|
||||
|
||||
acb_poly_div_series(res, A, B, len, prec);
|
||||
|
||||
arb_const_pi(c, prec);
|
||||
_acb_vec_scalar_mul_arb(res->coeffs, res->coeffs, res->length, c, prec);
|
||||
|
||||
acb_poly_clear(s);
|
||||
acb_poly_clear(u);
|
||||
acb_poly_clear(A);
|
||||
acb_poly_clear(B);
|
||||
acb_poly_clear(aa + 0);
|
||||
acb_poly_clear(aa + 1);
|
||||
acb_poly_clear(aa + 2);
|
||||
arb_clear(c);
|
||||
}
|
||||
|
||||
void
|
||||
acb_hypgeom_u_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, long prec)
|
||||
{
|
||||
if (acb_is_int(b))
|
||||
{
|
||||
acb_poly_t aa, bb, zz;
|
||||
|
||||
acb_poly_init(aa);
|
||||
acb_poly_init(bb);
|
||||
acb_poly_init(zz);
|
||||
|
||||
acb_poly_set_acb(aa, a);
|
||||
acb_poly_set_coeff_acb(bb, 0, b);
|
||||
acb_poly_set_coeff_si(bb, 1, 1);
|
||||
acb_poly_set_acb(zz, z);
|
||||
|
||||
acb_hypgeom_u_1f1_series(zz, aa, bb, zz, 1, prec);
|
||||
acb_poly_get_coeff_acb(res, zz, 0);
|
||||
|
||||
acb_poly_clear(aa);
|
||||
acb_poly_clear(bb);
|
||||
acb_poly_clear(zz);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_t t, u, v;
|
||||
acb_struct aa[3];
|
||||
|
||||
acb_init(t);
|
||||
acb_init(u);
|
||||
acb_init(v);
|
||||
acb_init(aa + 0);
|
||||
acb_init(aa + 1);
|
||||
acb_init(aa + 2);
|
||||
|
||||
acb_set(aa, a);
|
||||
acb_set(aa + 1, b);
|
||||
acb_one(aa + 2);
|
||||
acb_hypgeom_pfq_direct(u, aa, 1, aa + 1, 2, z, -1, prec);
|
||||
acb_sub(aa, a, b, prec);
|
||||
acb_add_ui(aa, aa, 1, prec);
|
||||
acb_sub_ui(aa + 1, b, 2, prec);
|
||||
acb_neg(aa + 1, aa + 1);
|
||||
acb_hypgeom_pfq_direct(v, aa, 1, aa + 1, 2, z, -1, prec);
|
||||
|
||||
acb_sub_ui(aa + 1, b, 1, prec);
|
||||
|
||||
/* rgamma(a-b+1) * gamma(1-b) * u */
|
||||
acb_rgamma(t, aa, prec);
|
||||
acb_mul(u, u, t, prec);
|
||||
acb_neg(t, aa + 1);
|
||||
acb_gamma(t, t, prec);
|
||||
acb_mul(u, u, t, prec);
|
||||
|
||||
/* rgamma(a) * gamma(b-1) * z^(1-b) * v */
|
||||
acb_rgamma(t, a, prec);
|
||||
acb_mul(v, v, t, prec);
|
||||
acb_gamma(t, aa + 1, prec);
|
||||
acb_mul(v, v, t, prec);
|
||||
acb_neg(t, aa + 1);
|
||||
acb_pow(t, z, t, prec);
|
||||
acb_mul(v, v, t, prec);
|
||||
|
||||
acb_add(res, u, v, prec);
|
||||
|
||||
acb_clear(t);
|
||||
acb_clear(u);
|
||||
acb_clear(v);
|
||||
acb_clear(aa + 0);
|
||||
acb_clear(aa + 1);
|
||||
acb_clear(aa + 2);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
acb_hypgeom_u(acb_t res, const acb_t a, const acb_t b, const acb_t z, long prec)
|
||||
{
|
||||
acb_t t;
|
||||
acb_init(t);
|
||||
|
||||
acb_sub(t, a, b, prec);
|
||||
acb_add_ui(t, t, 1, prec);
|
||||
|
||||
if ((acb_is_int(a) && arf_sgn(arb_midref(acb_realref(a))) <= 0) ||
|
||||
(acb_is_int(t) && arf_sgn(arb_midref(acb_realref(t))) <= 0) ||
|
||||
acb_hypgeom_u_use_asymp(z, prec))
|
||||
{
|
||||
acb_neg(t, a);
|
||||
acb_pow(t, z, t, prec);
|
||||
acb_hypgeom_u_asymp(res, a, b, z, -1, prec);
|
||||
acb_mul(res, res, t, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_hypgeom_u_1f1(res, a, b, z, prec);
|
||||
}
|
||||
|
||||
acb_clear(t);
|
||||
}
|
||||
|
||||
|
|
@ -223,6 +223,31 @@ in terms of the following auxiliary quantities
|
|||
to evaluate `U(a,b,z)` to *prec* accurate bits (assuming that *a* and *b*
|
||||
are small).
|
||||
|
||||
.. function:: void acb_hypgeom_u_1f1_series(acb_poly_t res, const acb_poly_t a, const acb_poly_t b, const acb_poly_t z, long len, long prec)
|
||||
|
||||
Computes `U(a,b,z)` as a power series truncated to length *len*,
|
||||
given `a, b, z \in \mathbb{C}[[x]]`.
|
||||
If `b \in \mathbb{Z}`, it computes one extra derivative and removes
|
||||
the singularity.
|
||||
As currently implemented, the output is indeterminate if `b` is nonexact
|
||||
and contains an integer.
|
||||
|
||||
.. function:: void acb_hypgeom_u_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, long prec)
|
||||
|
||||
Computes `U(a,b,z)` as a sum of two convergent hypergeometric series.
|
||||
If `b \in \mathbb{Z}`, it computes
|
||||
the limit value via :func:`acb_hypgeom_u_1f1_series`.
|
||||
As currently implemented, the output is indeterminate if `b` is nonexact
|
||||
and contains an integer.
|
||||
|
||||
.. function:: void acb_hypgeom_u(acb_t res, const acb_t a, const acb_t b, const acb_t z, long prec)
|
||||
|
||||
Computes `U(a,b,z)` using an automatic algorithm choice. The
|
||||
function :func:`acb_hypgeom_u_asymp` is used
|
||||
if `a` or `a-b+1` is a nonpositive integer (in which
|
||||
case the asymptotic series terminates), or if *z* is sufficiently large.
|
||||
Otherwise :func:`acb_hypgeom_u_1f1` is used.
|
||||
|
||||
The error function
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue