mirror of
https://github.com/vale981/arb
synced 2025-03-05 09:21:38 -05:00
MAINT: use acb_hypgeom_m for lower incomplete gamma
This commit is contained in:
alex
parent
0f280312a8
commit
53fb93bbdd
5 changed files with 111 additions and 163 deletions
|
|
@ -126,8 +126,6 @@ void acb_hypgeom_gamma_upper_1f1b(acb_t res, const acb_t s, const acb_t z, int m
|
|||
void acb_hypgeom_gamma_upper_singular(acb_t res, slong s, const acb_t z, int modified, slong prec);
|
||||
void acb_hypgeom_gamma_upper(acb_t res, const acb_t s, const acb_t z, int modified, slong prec);
|
||||
|
||||
void acb_hypgeom_gamma_lower_1f1a(acb_t res, const acb_t s, const acb_t z, int modified, slong prec);
|
||||
void acb_hypgeom_gamma_lower_1f1b(acb_t res, const acb_t s, const acb_t z, int modified, slong prec);
|
||||
void acb_hypgeom_gamma_lower(acb_t res, const acb_t s, const acb_t z, int modified, slong prec);
|
||||
|
||||
void acb_hypgeom_expint(acb_t res, const acb_t s, const acb_t z, slong prec);
|
||||
|
|
|
|||
|
|
@ -25,95 +25,46 @@
|
|||
|
||||
#include "acb_hypgeom.h"
|
||||
|
||||
void
|
||||
acb_hypgeom_gamma_lower_1f1a(acb_t res, const acb_t s,
|
||||
const acb_t z, int regularized, slong prec)
|
||||
{
|
||||
acb_t a, w;
|
||||
acb_struct b[2];
|
||||
|
||||
acb_init(a);
|
||||
acb_init(b);
|
||||
acb_init(b + 1);
|
||||
acb_init(w);
|
||||
|
||||
acb_set(a, s);
|
||||
acb_add_ui(b, s, 1, prec);
|
||||
acb_one(b + 1);
|
||||
acb_neg(w, z);
|
||||
|
||||
/* res = 1F1(s, s+1, -z) / s */
|
||||
acb_hypgeom_pfq_direct(res, a, 1, b, 2, w, -1, prec);
|
||||
acb_div(res, res, s, prec);
|
||||
|
||||
if (regularized == 0)
|
||||
{
|
||||
acb_pow(a, z, s, prec);
|
||||
acb_mul(res, res, a, prec);
|
||||
}
|
||||
else if (regularized == 1)
|
||||
{
|
||||
acb_pow(a, z, s, prec);
|
||||
acb_mul(res, res, a, prec);
|
||||
acb_rgamma(a, s, prec);
|
||||
acb_mul(res, res, a, prec);
|
||||
}
|
||||
|
||||
acb_clear(a);
|
||||
acb_clear(b);
|
||||
acb_clear(b + 1);
|
||||
acb_clear(w);
|
||||
}
|
||||
|
||||
void
|
||||
acb_hypgeom_gamma_lower_1f1b(acb_t res, const acb_t s,
|
||||
const acb_t z, int regularized, slong prec)
|
||||
{
|
||||
acb_t a, b;
|
||||
|
||||
acb_init(a);
|
||||
acb_init(b);
|
||||
|
||||
acb_add_ui(b, s, 1, prec);
|
||||
acb_hypgeom_pfq_direct(res, NULL, 0, b, 1, z, -1, prec);
|
||||
acb_div(res, res, s, prec);
|
||||
|
||||
acb_neg(a, z);
|
||||
acb_exp(a, a, prec);
|
||||
acb_mul(res, res, a, prec);
|
||||
|
||||
if (regularized == 0)
|
||||
{
|
||||
acb_pow(a, z, s, prec);
|
||||
acb_mul(res, res, a, prec);
|
||||
}
|
||||
else if (regularized == 1)
|
||||
{
|
||||
acb_pow(a, z, s, prec);
|
||||
acb_mul(res, res, a, prec);
|
||||
acb_rgamma(a, s, prec);
|
||||
acb_mul(res, res, a, prec);
|
||||
}
|
||||
|
||||
acb_clear(a);
|
||||
acb_clear(b);
|
||||
}
|
||||
|
||||
void
|
||||
acb_hypgeom_gamma_lower(acb_t res, const acb_t s, const acb_t z, int regularized, slong prec)
|
||||
{
|
||||
acb_t s1, nz, t;
|
||||
|
||||
if (!acb_is_finite(s) || !acb_is_finite(z))
|
||||
{
|
||||
acb_indeterminate(res);
|
||||
return;
|
||||
}
|
||||
|
||||
/* todo: handle the case where z is zero */
|
||||
acb_init(s1);
|
||||
acb_init(nz);
|
||||
acb_init(t);
|
||||
|
||||
/* todo: handle the case where s is an integer */
|
||||
acb_add_ui(s1, s, 1, prec);
|
||||
acb_neg(nz, z);
|
||||
|
||||
if (arf_sgn(arb_midref(acb_realref(z))) > 0)
|
||||
acb_hypgeom_gamma_lower_1f1b(res, s, z, regularized, prec);
|
||||
else
|
||||
acb_hypgeom_gamma_lower_1f1a(res, s, z, regularized, prec);
|
||||
if (regularized == 0)
|
||||
{
|
||||
/* \gamma(s, z) = s^-1 z^s 1F1(s, 1+s, -z) */
|
||||
acb_hypgeom_m(res, s, s1, nz, 0, prec);
|
||||
acb_pow(t, z, s, prec);
|
||||
acb_mul(res, res, t, prec);
|
||||
acb_div(res, res, s, prec);
|
||||
}
|
||||
else if (regularized == 1)
|
||||
{
|
||||
/* P(s, z) = z^s \gamma^{*}(s, z) */
|
||||
acb_hypgeom_m(res, s, s1, nz, 1, prec);
|
||||
acb_pow(t, z, s, prec);
|
||||
acb_mul(res, res, t, prec);
|
||||
}
|
||||
else if (regularized == 2)
|
||||
{
|
||||
/* \gamma^{*}(s, z) */
|
||||
acb_hypgeom_m(res, s, s1, nz, 1, prec);
|
||||
}
|
||||
|
||||
acb_clear(s1);
|
||||
acb_clear(nz);
|
||||
acb_clear(t);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -37,7 +37,7 @@ int main()
|
|||
|
||||
for (iter = 0; iter < 2000 * arb_test_multiplier(); iter++)
|
||||
{
|
||||
acb_t a0, a1, b, z, w0, w1, t, u;
|
||||
acb_t a0, a1, b, z, w0, w1, t, u, enz;
|
||||
slong prec0, prec1;
|
||||
int regularized;
|
||||
|
||||
|
|
@ -49,6 +49,7 @@ int main()
|
|||
acb_init(w1);
|
||||
acb_init(t);
|
||||
acb_init(u);
|
||||
acb_init(enz);
|
||||
|
||||
regularized = n_randint(state, 3);
|
||||
|
||||
|
|
@ -62,95 +63,45 @@ int main()
|
|||
|
||||
acb_add_ui(a1, a0, 1, prec0);
|
||||
|
||||
switch (n_randint(state, 3))
|
||||
{
|
||||
case 0:
|
||||
acb_hypgeom_gamma_lower_1f1a(w0, a0, z, regularized, prec0);
|
||||
break;
|
||||
case 1:
|
||||
acb_hypgeom_gamma_lower_1f1b(w0, a0, z, regularized, prec0);
|
||||
break;
|
||||
default:
|
||||
acb_hypgeom_gamma_lower(w0, a0, z, regularized, prec0);
|
||||
}
|
||||
acb_hypgeom_gamma_lower(w0, a0, z, regularized, prec0);
|
||||
acb_hypgeom_gamma_lower(w1, a1, z, regularized, prec1);
|
||||
|
||||
switch (n_randint(state, 3))
|
||||
{
|
||||
case 0:
|
||||
acb_hypgeom_gamma_lower_1f1a(w1, a0, z, regularized, prec1);
|
||||
break;
|
||||
case 1:
|
||||
acb_hypgeom_gamma_lower_1f1b(w1, a0, z, regularized, prec1);
|
||||
break;
|
||||
default:
|
||||
acb_hypgeom_gamma_lower(w1, a0, z, regularized, prec1);
|
||||
}
|
||||
|
||||
if (!acb_overlaps(w0, w1))
|
||||
{
|
||||
flint_printf("FAIL: consistency\n\n");
|
||||
flint_printf("a0 = "); acb_printd(a0, 30); flint_printf("\n\n");
|
||||
flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n");
|
||||
flint_printf("w0 = "); acb_printd(w0, 30); flint_printf("\n\n");
|
||||
flint_printf("w1 = "); acb_printd(w1, 30); flint_printf("\n\n");
|
||||
abort();
|
||||
}
|
||||
|
||||
switch (n_randint(state, 3))
|
||||
{
|
||||
case 1:
|
||||
acb_hypgeom_gamma_lower_1f1a(w1, a1, z, regularized, prec1);
|
||||
break;
|
||||
case 2:
|
||||
acb_hypgeom_gamma_lower_1f1b(w1, a1, z, regularized, prec1);
|
||||
break;
|
||||
default:
|
||||
acb_hypgeom_gamma_lower(w1, a1, z, regularized, prec1);
|
||||
}
|
||||
acb_neg(enz, z);
|
||||
acb_exp(enz, enz, prec0);
|
||||
|
||||
/* recurrence relations */
|
||||
if (regularized == 2)
|
||||
{
|
||||
/* a r(a,z) - exp(-z) - z r(a+1,z) = 0 */
|
||||
acb_neg(u, z);
|
||||
acb_exp(u, u, prec0);
|
||||
acb_neg(t, u);
|
||||
|
||||
acb_mul(b, w1, z, prec0);
|
||||
acb_addmul(t, a0, w0, prec0);
|
||||
acb_sub(t, t, b, prec0);
|
||||
/* gamma^{*}(a,z) - exp(-z)/Gamma(a+1) - z gamma^{*}(a+1,z) = 0 */
|
||||
/* http://dlmf.nist.gov/8.8.E4 */
|
||||
acb_set(t, w0);
|
||||
acb_rgamma(u, a1, prec0);
|
||||
acb_submul(t, enz, u, prec0);
|
||||
acb_submul(t, z, w1, prec0);
|
||||
}
|
||||
else if (regularized == 1)
|
||||
{
|
||||
/* q(a,z) - exp(-z) z^a / Gamma(a+1) - q(a+1,z) = 0 */
|
||||
acb_pow(t, z, a0, prec0);
|
||||
acb_rgamma(u, a1, prec0);
|
||||
acb_mul(t, t, u, prec0);
|
||||
|
||||
acb_neg(u, z);
|
||||
acb_exp(u, u, prec0);
|
||||
acb_neg(u, u);
|
||||
acb_mul(t, t, u, prec0);
|
||||
|
||||
acb_add(t, t, w0, prec0);
|
||||
acb_sub(t, t, w1, prec0);
|
||||
/* P(a,z) - exp(-z) z^a / Gamma(a+1) - P(a+1,z) = 0 */
|
||||
/* http://dlmf.nist.gov/8.8.E5 */
|
||||
acb_pow(u, z, a0, prec0);
|
||||
acb_rgamma(b, a1, prec0);
|
||||
acb_mul(u, u, b, prec0);
|
||||
acb_sub(t, w0, w1, prec0);
|
||||
acb_submul(t, enz, u, prec0);
|
||||
}
|
||||
else
|
||||
{
|
||||
/* a gamma(a,z) - exp(-z) z^a - gamma(a+1,z) = 0 */
|
||||
acb_pow(t, z, a0, prec0);
|
||||
|
||||
acb_neg(u, z);
|
||||
acb_exp(u, u, prec0);
|
||||
acb_neg(u, u);
|
||||
acb_mul(t, t, u, prec0);
|
||||
|
||||
acb_addmul(t, a0, w0, prec0);
|
||||
/* http://dlmf.nist.gov/8.8.E1 */
|
||||
acb_pow(u, z, a0, prec0);
|
||||
acb_mul(t, a0, w0, prec0);
|
||||
acb_submul(t, enz, u, prec0);
|
||||
acb_sub(t, t, w1, prec0);
|
||||
}
|
||||
|
||||
if (!acb_contains_zero(t))
|
||||
{
|
||||
flint_printf("FAIL: contiguous relation\n\n");
|
||||
flint_printf("FAIL: recurrence relation\n\n");
|
||||
flint_printf("regularized = %d\n\n", regularized);
|
||||
flint_printf("a0 = "); acb_printd(a0, 30); flint_printf("\n\n");
|
||||
flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n");
|
||||
|
|
@ -160,6 +111,50 @@ int main()
|
|||
abort();
|
||||
}
|
||||
|
||||
/* identities relating lower and upper incomplete gamma functions */
|
||||
if (regularized == 0 || regularized == 1)
|
||||
{
|
||||
acb_t u0;
|
||||
acb_init(u0);
|
||||
acb_hypgeom_gamma_upper(u0, a0, z, regularized, prec0);
|
||||
|
||||
acb_zero(t);
|
||||
|
||||
if (regularized == 1)
|
||||
{
|
||||
/* P(s,z) + Q(s,z) - 1 = 0 */
|
||||
/* http://dlmf.nist.gov/8.2.E5 */
|
||||
acb_add(t, w0, u0, prec0);
|
||||
acb_sub_ui(t, t, 1, prec0);
|
||||
}
|
||||
else
|
||||
{
|
||||
/* gamma(s,z) + Gamma(s,z) - Gamma(s) = 0 */
|
||||
/* excludes non-positive integer values of s */
|
||||
/* http://dlmf.nist.gov/8.2.E3 */
|
||||
if (!acb_is_int(a0) || arb_is_positive(acb_realref(a0)))
|
||||
{
|
||||
acb_gamma(b, a0, prec0);
|
||||
acb_add(t, w0, u0, prec0);
|
||||
acb_sub(t, t, b, prec0);
|
||||
}
|
||||
}
|
||||
|
||||
if (!acb_contains_zero(t))
|
||||
{
|
||||
flint_printf("FAIL: lower plus upper\n\n");
|
||||
flint_printf("regularized = %d\n\n", regularized);
|
||||
flint_printf("a0 = "); acb_printd(a0, 30); flint_printf("\n\n");
|
||||
flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n");
|
||||
flint_printf("w0 = "); acb_printd(w0, 30); flint_printf("\n\n");
|
||||
flint_printf("w1 = "); acb_printd(w1, 30); flint_printf("\n\n");
|
||||
flint_printf("t = "); acb_printd(t, 30); flint_printf("\n\n");
|
||||
abort();
|
||||
}
|
||||
|
||||
acb_clear(u0);
|
||||
}
|
||||
|
||||
acb_clear(a0);
|
||||
acb_clear(a1);
|
||||
acb_clear(b);
|
||||
|
|
@ -168,6 +163,7 @@ int main()
|
|||
acb_clear(w1);
|
||||
acb_clear(t);
|
||||
acb_clear(u);
|
||||
acb_clear(enz);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
|
|
|
|||
|
|
@ -177,6 +177,7 @@ int main()
|
|||
else if (regularized == 1)
|
||||
{
|
||||
/* Q(a,z) + exp(-z) z^a / Gamma(a+1) - Q(a+1,z) = 0 */
|
||||
/* http://dlmf.nist.gov/8.8.E6 */
|
||||
acb_pow(t, z, a0, prec0);
|
||||
acb_rgamma(u, a1, prec0);
|
||||
acb_mul(t, t, u, prec0);
|
||||
|
|
@ -191,6 +192,7 @@ int main()
|
|||
else
|
||||
{
|
||||
/* a Gamma(a,z) + exp(-z) z^a - Gamma(a+1,z) = 0 */
|
||||
/* http://dlmf.nist.gov/8.8.E2 */
|
||||
acb_pow(t, z, a0, prec0);
|
||||
|
||||
acb_neg(u, z);
|
||||
|
|
|
|||
|
|
@ -589,15 +589,16 @@ 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_lower_1f1a(acb_t res, const acb_t s, const acb_t z, int regularized, slong prec)
|
||||
|
||||
.. function:: void acb_hypgeom_gamma_lower_1f1b(acb_t res, const acb_t s, const acb_t z, int regularized, slong prec)
|
||||
|
||||
.. function:: void acb_hypgeom_gamma_lower(acb_t res, const acb_t s, const acb_t z, int regularized, slong prec)
|
||||
|
||||
Computes the lower incomplete gamma function `\gamma(s,z)`
|
||||
with regularization options similar to those for the upper
|
||||
incomplete gamma function.
|
||||
If *regularized* is 0, computes the lower incomplete gamma function
|
||||
`\gamma(s,z)`.
|
||||
|
||||
If *regularized* is 1, computes the regularized lower incomplete
|
||||
gamma function `P(s,z) = \Gamma(s,z) / \Gamma(s)`.
|
||||
|
||||
If *regularized* is 2, computes a further regularized lower incomplete
|
||||
gamma function `\gamma^{*}(s,z) = z^{-s} P(s,z)`.
|
||||
|
||||
Exponential and trigonometric integrals
|
||||
-------------------------------------------------------------------------------
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue