Skip to content

Commit

Permalink
Merge pull request flintlib#2130 from fredrik-johansson/eintegrals
Browse files Browse the repository at this point in the history 
Good enclosure for exponential and trigonometric integrals
  • Loading branch information
@fredrik-johansson
fredrik-johansson authored Dec 30, 2024
2 parents 972b3fa + 74014f1 commit c0a82ba
Show file tree
Hide file tree
Showing 6 changed files with 469 additions and 64 deletions.
21 changes: 20 additions & 1 deletion examples/integrals.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -732,11 +732,22 @@ scaled_bessel_select_N(arb_t N, ulong k, slong prec)
arb_mul_2exp_si(N, N, e);
}

/* f(z) = Si(z) */
int
f_si(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */

acb_hypgeom_si(res, z, prec);
return 0;
}

/* ------------------------------------------------------------------------- */
/* Main test program */
/* ------------------------------------------------------------------------- */

#define NUM_INTEGRALS 37
#define NUM_INTEGRALS 38

const char * descr[NUM_INTEGRALS] =
{
Expand Down Expand Up @@ -777,6 +788,7 @@ const char * descr[NUM_INTEGRALS] =
"int_0^{inf} exp(-x) I_0(x/15)^{15} dx (using domain truncation)",
"int_{-1-i}^{-1+i} 1/sqrt(x) dx",
"int_0^{inf} 1/gamma(x) dx (using domain truncation)",
"int_0^{1000} Si(x) dx",
};

int main(int argc, char *argv[])
Expand Down Expand Up @@ -1273,6 +1285,13 @@ int main(int argc, char *argv[])
arb_add_error_2exp_si(acb_realref(s), -goal);
break;

case 37:
acb_zero(a);
acb_set_ui(b, 1000);
acb_calc_integrate(s, f_si, NULL, a, b, goal, tol, options, prec);
break;


default:
abort();
}
Expand Down
141 changes: 105 additions & 36 deletions src/acb_hypgeom/chi.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
(at your option) any later version. See <https://www.gnu.org/licenses/>.
*/

#include <math.h>
#include "acb.h"
#include "acb_hypgeom.h"

Expand Down Expand Up @@ -101,51 +102,119 @@ acb_hypgeom_chi_asymp(acb_t res, const acb_t z, slong prec)
acb_clear(one);
}

/* defined in ci.c */
void
acb_hypgeom_chi_2f3(acb_t res, const acb_t z, slong prec)
{
acb_t a, t, u;
acb_struct b[3];

acb_init(a);
acb_init(b);
acb_init(b + 1);
acb_init(b + 2);
acb_init(t);
acb_init(u);

acb_one(a);
acb_set_ui(b, 2);
acb_set(b + 1, b);
acb_set_ui(b + 2, 3);
acb_mul_2exp_si(b + 2, b + 2, -1);

acb_mul(t, z, z, prec);
acb_mul_2exp_si(t, t, -2);
acb_hypgeom_pfq_direct(u, a, 1, b, 3, t, -1, prec);
acb_mul(u, u, t, prec);
acb_hypgeom_chi_2f3(acb_t res, const acb_t z, slong prec);

acb_log(t, z, prec);
acb_add(u, u, t, prec);
/* Bound propagated error |chi(z)-chi(mid(z))| assuming that we
are off the branch cut (not checked here).
Uses |chi'(z)| = |cosh(z)/z| <= cosh(|re(z)|)/|z|.
*/
static void
acb_hypgeom_chi_prop_error(mag_t err, const acb_t z)
{
if (acb_is_exact(z))
{
mag_zero(err);
}
else
{
mag_t t, u;

arb_const_euler(acb_realref(t), prec);
arb_add(acb_realref(u), acb_realref(u), acb_realref(t), prec);
mag_init(t);
mag_init(u);

acb_swap(res, u);
arb_get_mag(t, acb_realref(z));
mag_cosh(t, t);
acb_get_mag_lower(u, z);
mag_div(t, t, u);
mag_hypot(u, arb_radref(acb_realref(z)), arb_radref(acb_imagref(z)));
mag_mul(err, t, u);

acb_clear(a);
acb_clear(b);
acb_clear(b + 1);
acb_clear(b + 2);
acb_clear(t);
acb_clear(u);
mag_clear(t);
mag_clear(u);
}
}

void
acb_hypgeom_chi(acb_t res, const acb_t z, slong prec)
{
if (acb_hypgeom_u_use_asymp(z, prec))
acb_hypgeom_chi_asymp(res, z, prec);
if (!acb_is_finite(z) || acb_contains_zero(z))
{
acb_indeterminate(res);
}
/* todo -- good real code
else if (acb_is_real(z))
{
}
*/
else if (arb_contains_zero(acb_imagref(z)) &&
!arb_is_nonnegative(acb_imagref(z)) &&
!arb_is_positive(acb_realref(z)))
{
/* We straddle the branch cut; do something generic */
if (acb_hypgeom_u_use_asymp(z, prec))
acb_hypgeom_chi_asymp(res, z, prec);
else
acb_hypgeom_chi_2f3(res, z, prec);
}
else
acb_hypgeom_chi_2f3(res, z, prec);
{
double asymp_accuracy, a, b, absz, cancellation, rlog2 = 1.4426950408889634;
slong wp;
int use_asymp;
acb_t m;
mag_t err;
/* since we don't have a special case for real z */
int pure_real = arb_is_zero(acb_imagref(z));

a = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
b = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_DOWN);
a = fabs(a);
b = fabs(b);
absz = sqrt(a * a + b * b);

if (a <= 1.0 && b <= 1.0)
{
use_asymp = 0;
}
else if (a > prec || b > prec)
{
use_asymp = 1;
}
else
{
asymp_accuracy = absz * rlog2 * 0.999 - 5;
use_asymp = asymp_accuracy > prec;
}

acb_init(m);
mag_init(err);

/* assumes that we already handled the branch cut */
acb_hypgeom_chi_prop_error(err, z);
acb_get_mid(m, z);

if (use_asymp)
{
acb_hypgeom_chi_asymp(res, m, prec);
}
else
{
/* terms grow to ~ exp(|z|), sum is ~ exp(|re(z)|) */
cancellation = (absz - a) * rlog2;
wp = prec + FLINT_MAX(0, cancellation);
wp = wp * 1.001 + 5;
acb_hypgeom_chi_2f3(res, m, wp);
acb_set_round(res, res, prec);
}

if (pure_real)
arb_add_error_mag(acb_realref(res), err);
else
acb_add_error_mag(res, err);

acb_clear(m);
mag_clear(err);
}
}
128 changes: 115 additions & 13 deletions src/acb_hypgeom/ci.c
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
/*
Copyright (C) 2015 Fredrik Johansson
Copyright (C) 2015, 2024 Fredrik Johansson
This file is part of FLINT.
Expand All @@ -9,6 +9,7 @@
(at your option) any later version. See <https://www.gnu.org/licenses/>.
*/

#include <math.h>
#include "acb.h"
#include "arb_hypgeom.h"
#include "acb_hypgeom.h"
Expand Down Expand Up @@ -109,8 +110,8 @@ acb_hypgeom_ci_asymp(acb_t res, const acb_t z, slong prec)
acb_clear(one);
}

void
acb_hypgeom_ci_2f3(acb_t res, const acb_t z, slong prec)
static void
_acb_hypgeom_ci_2f3(acb_t res, const acb_t z, int hyperbolic, slong prec)
{
acb_t a, t, u;
acb_struct b[3];
Expand All @@ -130,7 +131,8 @@ acb_hypgeom_ci_2f3(acb_t res, const acb_t z, slong prec)

acb_mul(t, z, z, prec);
acb_mul_2exp_si(t, t, -2);
acb_neg(t, t);
if (!hyperbolic)
acb_neg(t, t);
acb_hypgeom_pfq_direct(u, a, 1, b, 3, t, -1, prec);
acb_mul(u, u, t, prec);

Expand All @@ -150,31 +152,131 @@ acb_hypgeom_ci_2f3(acb_t res, const acb_t z, slong prec)
acb_clear(u);
}

void
acb_hypgeom_ci_2f3(acb_t res, const acb_t z, slong prec)
{
_acb_hypgeom_ci_2f3(res, z, 0, prec);
}

void
acb_hypgeom_chi_2f3(acb_t res, const acb_t z, slong prec)
{
_acb_hypgeom_ci_2f3(res, z, 1, prec);
}

/* Bound propagated error |ci(z)-ci(mid(z))| assuming that we
are off the branch cut (not checked here).
Uses |ci'(z)| = |cos(z)/z| <= cosh(|im(z)|)/|z|.
*/
static void
acb_hypgeom_ci_prop_error(mag_t err, const acb_t z)
{
if (acb_is_exact(z))
{
mag_zero(err);
}
else
{
mag_t t, u;

mag_init(t);
mag_init(u);

arb_get_mag(t, acb_imagref(z));
mag_cosh(t, t);
acb_get_mag_lower(u, z);
mag_div(t, t, u);
mag_hypot(u, arb_radref(acb_realref(z)), arb_radref(acb_imagref(z)));
mag_mul(err, t, u);

mag_clear(t);
mag_clear(u);
}
}

void
acb_hypgeom_ci(acb_t res, const acb_t z, slong prec)
{
if (acb_is_real(z) && arb_is_finite(acb_realref(z)))
if (!acb_is_finite(z) || acb_contains_zero(z))
{
acb_indeterminate(res);
}
else if (acb_is_real(z))
{
if (arb_is_positive(acb_realref(z)))
{
arb_hypgeom_ci(acb_realref(res), acb_realref(z), prec);
arb_zero(acb_imagref(res));
}
else if (arb_is_negative(acb_realref(z)))
else /* arb_is_negative(acb_realref(z)), since we already tested for 0 */
{
arb_neg(acb_realref(res), acb_realref(z));
arb_hypgeom_ci(acb_realref(res), acb_realref(res), prec);
arb_const_pi(acb_imagref(res), prec);
}
}
else if (arb_contains_zero(acb_imagref(z)) &&
!arb_is_nonnegative(acb_imagref(z)) &&
!arb_is_positive(acb_realref(z)))
{
/* We straddle the branch cut; do something generic */
if (acb_hypgeom_u_use_asymp(z, prec))
acb_hypgeom_ci_asymp(res, z, prec);
else
acb_hypgeom_ci_2f3(res, z, prec);
}
else
{
double asymp_accuracy, a, b, absz, cancellation, rlog2 = 1.4426950408889634;
slong wp;
int use_asymp;
acb_t m;
mag_t err;

a = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
b = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_DOWN);
a = fabs(a);
b = fabs(b);
absz = sqrt(a * a + b * b);

if (a <= 1.0 && b <= 1.0)
{
acb_indeterminate(res);
use_asymp = 0;
}
else if (a > prec || b > prec)
{
use_asymp = 1;
}
else
{
asymp_accuracy = absz * rlog2 * 0.999 - 5;
use_asymp = asymp_accuracy > prec;
}
return;
}

if (acb_hypgeom_u_use_asymp(z, prec))
acb_hypgeom_ci_asymp(res, z, prec);
else
acb_hypgeom_ci_2f3(res, z, prec);
acb_init(m);
mag_init(err);

/* assumes that we already handled the branch cut */
acb_hypgeom_ci_prop_error(err, z);
acb_get_mid(m, z);

if (use_asymp)
{
acb_hypgeom_ci_asymp(res, m, prec);
}
else
{
/* terms grow to ~ exp(|z|), sum is ~ exp(|im(z)|) */
cancellation = (absz - b) * rlog2;
wp = prec + FLINT_MAX(0, cancellation);
wp = wp * 1.001 + 5;
acb_hypgeom_ci_2f3(res, m, wp);
acb_set_round(res, res, prec);
}

acb_add_error_mag(res, err);

acb_clear(m);
mag_clear(err);
}
}
Loading

0 comments on commit c0a82ba

Please sign in to comment.