hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-04 17:01:40 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions
arb/acb_hypgeom/u.c

473 lines
12 KiB
C
Raw Permalink Blame History

/*
Copyright (C) 2015 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 "arb_hypgeom.h"
#include "acb_hypgeom.h"
static void
bsplit(acb_t A, acb_t B, acb_t C, acb_t D,
const acb_t b, const acb_t z, slong n0, slong n1, slong prec)
{
if (n1 - n0 == 1)
{
acb_zero(A);
acb_one(B);
acb_neg(C, b);
acb_add_si(C, C, 2 - n0, prec);
acb_mul_si(C, C, n0 - 1, prec);
acb_sub(D, z, b, prec);
acb_add_si(D, D, 2 - 2 * n0, prec);
}
else
{
slong m;
acb_t T, A2, B2, C2, D2;
acb_init(T);
acb_init(A2);
acb_init(B2);
acb_init(C2);
acb_init(D2);
m = n0 + (n1 - n0) / 2;
bsplit(A, B, C, D, b, z, n0, m, prec);
bsplit(A2, B2, C2, D2, b, z, m, n1, prec);
acb_set(T, A);
acb_mul(A, A, A2, prec);
acb_addmul(A, B2, C, prec);
acb_mul(C, C, D2, prec);
acb_addmul(C, C2, T, prec);
acb_set(T, B);
acb_mul(B, B, A2, prec);
acb_addmul(B, B2, D, prec);
acb_mul(D, D, D2, prec);
acb_addmul(D, C2, T, prec);
acb_clear(T);
acb_clear(A2);
acb_clear(B2);
acb_clear(C2);
acb_clear(D2);
}
}
void
acb_hypgeom_u_si_rec(acb_t res, slong a, const acb_t b, const acb_t z, slong prec)
{
slong k;
acb_t u0, u1, t;
if (a > 0)
flint_abort();
if (a == 0)
{
acb_one(res);
return;
}
else if (a == -1)
{
acb_sub(res, z, b, prec);
return;
}
/* special-case U(-n, -n+1, z) = z^n */
if (acb_equal_si(b, a + 1))
{
acb_pow_si(res, z, -a, prec);
return;
}
acb_init(u0);
acb_init(u1);
acb_init(t);
acb_one(u0);
acb_sub(u1, z, b, prec);
if (-a < 20)
{
for (k = 2; k <= -a; k++)
{
acb_neg(t, b);
acb_add_si(t, t, 2 - k, prec);
acb_mul_si(t, t, k - 1, prec);
acb_mul(u0, u0, t, prec);
acb_sub(t, z, b, prec);
acb_add_si(t, t, 2 - 2 * k, prec);
acb_addmul(u0, u1, t, prec);
acb_swap(u0, u1);
}
acb_set(res, u1);
}
else
{
acb_t A, B, C, D;
acb_init(A);
acb_init(B);
acb_init(C);
acb_init(D);
bsplit(A, B, C, D, b, z, 2, -a + 1, prec);
acb_sub(A, z, b, prec);
acb_mul(D, D, A, prec);
acb_add(res, C, D, prec);
acb_clear(A);
acb_clear(B);
acb_clear(C);
acb_clear(D);
}
acb_clear(u0);
acb_clear(u1);
acb_clear(t);
}
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,
slong len, slong prec)
{
acb_poly_t s, u, A, B;
acb_poly_struct aa[3];
arb_t c;
slong 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, -1, wlen, prec);
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, -1, wlen, prec);
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, slong 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_choose(int * asymp, slong * wp,
const acb_t a, const acb_t b, const acb_t z, slong prec)
{
double x, y, t, cancellation;
double input_accuracy, direct_accuracy, asymp_accuracy;
*asymp = 0;
*wp = prec;
input_accuracy = acb_rel_one_accuracy_bits(z);
t = acb_rel_one_accuracy_bits(a);
input_accuracy = FLINT_MIN(input_accuracy, t);
t = acb_rel_one_accuracy_bits(b);
input_accuracy = FLINT_MIN(input_accuracy, t);
input_accuracy = FLINT_MAX(input_accuracy, 0.0);
/* From here we ignore the values of a, b. Taking them into account is
a possible future improvement... */
/* Tiny |z|. */
if ((arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 2) < 0 &&
arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), 2) < 0))
{
*asymp = 0;
*wp = FLINT_MAX(2, FLINT_MIN(input_accuracy + 20, prec));
return;
}
/* Huge |z|. */
if ((arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 64) > 0 ||
arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), 64) > 0))
{
*asymp = 1;
*wp = FLINT_MAX(2, FLINT_MIN(input_accuracy + 20, prec));
return;
}
x = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
y = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_DOWN);
asymp_accuracy = sqrt(x * x + y * y) * 1.44269504088896;
/* The Kummer transformation gives less cancellation with the 1F1 series.
if (x < 0.0)
{
*kummer = 1;
x = -x;
} */
if (asymp_accuracy >= prec)
{
*asymp = 1;
*wp = FLINT_MAX(2, FLINT_MIN(input_accuracy + 20, prec));
return;
}
/* Assume U ~ 1, M ~ exp(|z|) (there is cancellation both in the
evaluation of M and in the linear combination) -- a better estimate
would account for a, b. */
cancellation = sqrt(x * x + y * y) * 1.44269504088896 + 5;
direct_accuracy = input_accuracy - cancellation;
if (direct_accuracy > asymp_accuracy)
{
*asymp = 0;
*wp = FLINT_MAX(2, FLINT_MIN(input_accuracy + 20, prec + cancellation));
}
else
{
*asymp = 1;
*wp = FLINT_MAX(2, FLINT_MIN(input_accuracy + 20, prec));
}
}
void
acb_hypgeom_u_nointegration(acb_t res, const acb_t a, const acb_t b, const acb_t z, slong prec)
{
acb_t t;
arf_srcptr av, tv;
av = arb_midref(acb_realref(a));
/* Handle small polynomial cases without divisions. */
/* todo: should incorporate a -> 1+a-b transformation, also... */
if (acb_is_int(a) && arf_sgn(av) <= 0)
{
if (arf_cmpabs_ui(av, 30) < 0 ||
(arf_cmpabs_ui(av, prec) < 0 &&
((acb_bits(b) < 0.1 * prec && acb_bits(z) < 0.1 * prec)
|| acb_contains_zero(z))))
{
acb_hypgeom_u_si_rec(res, arf_get_si(av, ARF_RND_DOWN), b, z, prec);
return;
}
}
acb_init(t);
acb_sub(t, a, b, prec);
acb_add_ui(t, t, 1, prec);
tv = arb_midref(acb_realref(t));
/* todo: combine these conditions with the code below */
if ((acb_is_int(a) && arf_sgn(av) <= 0) ||
(acb_is_int(t) && arf_sgn(tv) <= 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
{
slong wp;
int asymp;
acb_hypgeom_u_choose(&asymp, &wp, a, b, z, prec);
if (asymp)
{
acb_neg(t, a);
acb_pow(t, z, t, prec);
acb_hypgeom_u_asymp(res, a, b, z, -1, wp);
acb_mul(res, res, t, prec);
}
else
{
acb_hypgeom_u_1f1(res, a, b, z, wp);
acb_set_round(res, res, prec);
}
}
acb_clear(t);
}
void
acb_hypgeom_u(acb_t res, const acb_t a, const acb_t b, const acb_t z, slong prec)
{
acb_t res2;
slong acc, max, t;
acb_init(res2);
acb_hypgeom_u_nointegration(res2, a, b, z, prec);
acc = acb_rel_accuracy_bits(res2);
if (acc < 0.5 * prec)
{
max = prec;
t = acb_rel_accuracy_bits(z);
max = FLINT_MIN(max, t);
t = acb_rel_accuracy_bits(a);
max = FLINT_MIN(max, t);
t = acb_rel_accuracy_bits(b);
max = FLINT_MIN(max, t);
if (max > 2 && acc < 0.5 * max)
{
if (acb_is_real(a) && acb_is_real(b) && acb_is_real(z) &&
arb_is_positive(acb_realref(z)) &&
arf_cmpabs_2exp_si(arb_midref(acb_realref(a)), 60) < 0 &&
arf_cmpabs_2exp_si(arb_midref(acb_realref(b)), 60) < 0 &&
arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 60) < 0)
{
arb_hypgeom_u_integration(acb_realref(res),
acb_realref(a), acb_realref(b), acb_realref(z), prec);
arb_zero(acb_imagref(res));
if (acb_rel_accuracy_bits(res) > acb_rel_accuracy_bits(res2) ||
(acb_is_finite(res) && !acb_is_finite(res2)))
{
acb_swap(res, res2);
}
}
}
}
acb_swap(res, res2);
acb_clear(res2);
}
Reference in a new issue View git blame Copy permalink