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

new arb_hypgeom_ci

Browse source
This commit is contained in:
fredrik 2021-12-21 22:53:06 +01:00
parent fcffcf32b8
commit 88ee2a68fb
8 changed files with 363 additions and 29 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

22
acb_hypgeom/ci.c
View file

@ -9,6 +9,7 @@
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "arb_hypgeom.h"
#include "acb_hypgeom.h"
void
@ -151,9 +152,28 @@ acb_hypgeom_ci_2f3(acb_t res, const acb_t z, slong prec)
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 (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)))
{
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
{
acb_indeterminate(res);
}
return;
}
if (acb_hypgeom_u_use_asymp(z, prec))
acb_hypgeom_ci_asymp(res, z, prec);
else
acb_hypgeom_ci_2f3(res, z, prec);
}

4
arb_hypgeom.h
View file

@ -100,7 +100,6 @@ void arb_hypgeom_ei(arb_t res, const arb_t z, slong prec);
void _arb_hypgeom_ei_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec);
void arb_hypgeom_ei_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec);
void _arb_hypgeom_si_asymp(arb_t res, const arb_t z, slong N, slong prec);
void _arb_hypgeom_si_1f2(arb_t res, const arb_t z, slong N, slong wp, slong prec);
void arb_hypgeom_si(arb_t res, const arb_t z, slong prec);
@ -108,7 +107,10 @@ void arb_hypgeom_si(arb_t res, const arb_t z, slong prec);
void _arb_hypgeom_si_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec);
void arb_hypgeom_si_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec);
void _arb_hypgeom_ci_asymp(arb_t res, const arb_t z, slong N, slong prec);
void _arb_hypgeom_ci_2f3(arb_t res, const arb_t z, slong N, slong wp, slong prec);
void arb_hypgeom_ci(arb_t res, const arb_t z, slong prec);
void _arb_hypgeom_ci_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec);
void arb_hypgeom_ci_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec);

230
arb_hypgeom/ci.c Normal file
View file

@ -0,0 +1,230 @@
/*
Copyright (C) 2021 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"
#define LOG2 0.69314718055994531
#define INV_LOG2 1.4426950408889634074
#define EXP1 2.7182818284590452
double arf_get_d_log2_abs_approx_clamped(const arf_t x);
void
_arb_hypgeom_ci_asymp(arb_t res, const arb_t z, slong N, slong prec)
{
arb_t s, c, sz, cz, u;
fmpq a[1];
slong wp;
mag_t err, t;
N = FLINT_MAX(N, 1);
arb_init(s);
arb_init(c);
arb_init(sz);
arb_init(cz);
arb_init(u);
mag_init(err);
mag_init(t);
/* Error is bounded by first omitted term, N! / z^N */
arb_get_mag_lower(err, z);
mag_inv(err, err);
mag_pow_ui(err, err, N);
mag_fac_ui(t, N);
mag_mul(err, err, t);
wp = prec * 1.001 + 5;
arb_set(u, z);
*fmpq_numref(&a[0]) = 1;
*fmpq_denref(&a[0]) = 1;
arb_hypgeom_sum_fmpq_imag_arb(c, s, a, 1, NULL, 0, u, 1, N, wp);
arb_add_error_mag(c, err);
arb_add_error_mag(s, err);
arb_sin_cos(sz, cz, z, wp);
arb_mul(c, c, sz, wp);
arb_submul(c, s, cz, wp);
arb_div(res, c, z, prec);
arb_clear(s);
arb_clear(c);
arb_clear(sz);
arb_clear(cz);
arb_clear(u);
mag_clear(err);
mag_clear(t);
}
void
_arb_hypgeom_ci_2f3(arb_t res, const arb_t z, slong N, slong wp, slong prec)
{
mag_t err, t;
arb_t s, u;
fmpq a[1];
fmpq b[3];
N = FLINT_MAX(N, 1);
mag_init(err);
mag_init(t);
arb_init(s);
arb_init(u);
arb_sqr(u, z, wp);
arb_mul_2exp_si(u, u, -2);
arb_neg(u, u);
*fmpq_numref(&a[0]) = 1;
*fmpq_denref(&a[0]) = 1;
*fmpq_numref(&b[0]) = 2;
*fmpq_denref(&b[0]) = 1;
*fmpq_numref(&b[1]) = 2;
*fmpq_denref(&b[1]) = 1;
*fmpq_numref(&b[2]) = 3;
*fmpq_denref(&b[2]) = 2;
/* Terms are bounded by u^N / (4 (N!)^2) */
arb_get_mag(err, u);
/* u^N */
mag_set(t, err);
mag_pow_ui(t, t, N);
/* geometric factor for u/N^2 */
mag_div_ui(err, err, N);
mag_div_ui(err, err, N);
mag_geom_series(err, err, 0);
mag_mul(t, t, err);
/* 1/(N!)^2 */
mag_rfac_ui(err, N);
mag_mul(err, err, err);
mag_mul(err, err, t);
/* 1/4 */
mag_mul_2exp_si(err, err, -2);
arb_hypgeom_sum_fmpq_arb(s, a, 1, b, 3, u, 0, N, wp);
arb_add_error_mag(s, err);
arb_mul(s, s, u, wp);
arb_log(u, z, wp);
arb_add(s, s, u, wp);
arb_const_euler(u, wp);
arb_add(res, s, u, prec);
mag_clear(err);
mag_clear(t);
arb_clear(u);
arb_clear(s);
}
void
arb_hypgeom_ci(arb_t res, const arb_t z, slong prec)
{
slong wp, N, acc;
double dz, du;
if (!arb_is_positive(z) || !arb_is_finite(z))
{
arb_indeterminate(res);
return;
}
if (ARF_IS_LAGOM(arb_midref(z)))
{
acc = arb_rel_accuracy_bits(z);
acc = FLINT_MAX(acc, 0);
acc = FLINT_MIN(acc, prec);
acc += FLINT_MAX(0, -ARF_EXP(arb_midref(z)));
prec = FLINT_MIN(prec, acc + 32);
}
dz = fabs(arf_get_d(arb_midref(z), ARF_RND_DOWN));
dz = FLINT_MIN(dz, 1e300);
if (dz > 2.0)
{
double log2_err, err_prev, log2dz;
log2dz = arf_get_d_log2_abs_approx_clamped(arb_midref(z));
err_prev = 0.0;
for (N = 1; N < 2 * prec; N++)
{
log2_err = ((N + 1.0) * (log(N + 1.0) - 1.0)) * INV_LOG2 - N * log2dz;
if (log2_err > err_prev)
break;
if (log2_err < -prec - 2)
{
_arb_hypgeom_ci_asymp(res, z, N, prec);
return;
}
err_prev = log2_err;
}
}
if (arf_cmpabs_2exp_si(arb_midref(z), -30) < 0)
{
N = -arf_abs_bound_lt_2exp_si(arb_midref(z));
wp = prec * 1.001 + 10;
N = (prec + N - 1) / N;
}
else
{
du = 0.25 * dz * dz;
wp = prec * 1.001 + 10;
if (du > 1.0)
wp += dz * 1.4426950408889634;
N = (prec + 5) * LOG2 / (2 * d_lambertw((prec + 5) * LOG2 / (2 * EXP1 * sqrt(du)))) + 1;
}
if (arb_is_exact(z))
{
_arb_hypgeom_ci_2f3(res, z, N, wp, prec);
}
else
{
mag_t err;
mag_init(err);
/* |ci'(z)| = |cos(z)/z| <= 1/z */
arb_get_mag_lower(err, z);
if (mag_cmp_2exp_si(err, 0) >= 0 || 1)
{
arb_t zmid;
arb_init(zmid);
arb_get_mid_arb(zmid, z);
mag_inv(err, err);
mag_mul(err, err, arb_radref(z));
_arb_hypgeom_ci_2f3(res, zmid, N, wp, prec);
arb_add_error_mag(res, err);
arb_clear(zmid);
}
else
{
_arb_hypgeom_ci_2f3(res, z, N, wp, prec);
}
mag_clear(err);
}
}

19
arb_hypgeom/sum_fmpq_arb_rs.c
View file

@ -34,12 +34,12 @@ tail_precision(slong k, double min_k, slong alen, slong blen, double log2z, doub
new_prec = FLINT_MIN(new_prec, prec);
new_prec = FLINT_MAX(new_prec, 32);
/* printf("term %ld, max %f, log2x %f, magn %f new_prec %ld\n", k, log2z, log2max, term_magnitude, new_prec); */
/* printf("term %ld, max %f, log2x %f, magn %f new_prec %ld\n", k, log2z, log2max, term_magnitude, new_prec); */
return new_prec;
}
/* Return approximation of log2(|x|), clambed between COEFF_MIN and COEFF_MAX. */
/* Return approximation of log2(|x|), approximately clamped between COEFF_MIN and COEFF_MAX. */
double
arf_get_d_log2_abs_approx_clamped(const arf_t x)
{
@ -60,12 +60,21 @@ arf_get_d_log2_abs_approx_clamped(const arf_t x)
}
else
{
mp_srcptr tp;
mp_size_t tn;
double v;
slong e = ARF_EXP(x);
if (e < -50 || e > 50)
return e;
ARF_GET_MPN_READONLY(tp, tn, x);
if (tn == 1)
v = (double)(tp[0]);
else
return 1.4426950408889634074 * mag_d_log_upper_bound(fabs(arf_get_d(x, ARF_RND_UP)));
v = (double)(tp[tn - 1]) + (double)(tp[tn - 2]) * ldexp(1.0, -FLINT_BITS);
v *= ldexp(1.0, -FLINT_BITS);
return 1.4426950408889634074 * mag_d_log_upper_bound(v) + (double) e;
}
}

92
arb_hypgeom/test/t-ci.c Normal file
View file

@ -0,0 +1,92 @@
/*
Copyright (C) 2018 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"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("ci....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
{
arb_t x, s, t;
slong prec1, prec2;
arb_init(x);
arb_init(s);
arb_init(t);
if (n_randint(state, 10) == 0)
{
prec1 = 2 + n_randint(state, 2000);
prec2 = 2 + n_randint(state, 2000);
}
else
{
prec1 = 2 + n_randint(state, 200);
prec2 = 2 + n_randint(state, 200);
}
arb_randtest(x, state, 2 + n_randint(state, prec1), 2 + n_randint(state, 100));
arb_randtest(s, state, 2 + n_randint(state, prec1), 2 + n_randint(state, 100));
arb_randtest(t, state, 2 + n_randint(state, prec1), 2 + n_randint(state, 100));
switch (n_randint(state, 3))
{
case 0:
arb_hypgeom_ci(s, x, prec1);
break;
case 1:
_arb_hypgeom_ci_2f3(s, x, n_randint(state, prec1), prec1, prec1);
break;
default:
_arb_hypgeom_c_asymp(s, x, n_randint(state, prec1 / 2), prec1);
break;
}
switch (n_randint(state, 3))
{
case 0:
arb_hypgeom_ci(t, x, prec2);
break;
case 1:
_arb_hypgeom_ci_2f3(t, x, n_randint(state, prec2), prec2, prec2);
break;
default:
_arb_hypgeom_ci_asymp(t, x, n_randint(state, prec2 / 2), prec2);
break;
}
if (!arb_overlaps(s, t))
{
flint_printf("FAIL: overlap\n\n");
flint_printf("x = "); arb_printn(x, 100, 0); flint_printf("\n\n");
flint_printf("s = "); arb_printn(s, 100, 0); flint_printf("\n\n");
flint_printf("t = "); arb_printn(t, 100, 0); flint_printf("\n\n");
flint_abort();
}
arb_clear(x);
arb_clear(s);
arb_clear(t);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

3
arb_hypgeom/test/t-si.c
View file

@ -71,9 +71,6 @@ int main()
break;
}
arb_hypgeom_si(s, x, prec1);
_arb_hypgeom_si_asymp(t, x, n_randint(state, prec2 / 2), prec2);
if (!arb_overlaps(s, t))
{
flint_printf("FAIL: overlap\n\n");

18
arb_hypgeom/wrappers.c
View file

@ -70,24 +70,6 @@ arb_hypgeom_ei(arb_t res, const arb_t z, slong prec)
}
}
void
arb_hypgeom_ci(arb_t res, const arb_t z, slong prec)
{
if (!arb_is_finite(z) || !arb_is_positive(z))
{
arb_indeterminate(res);
}
else
{
acb_t t;
acb_init(t);
arb_set(acb_realref(t), z);
acb_hypgeom_ci(t, t, prec);
arb_swap(res, acb_realref(t));
acb_clear(t);
}
}
void
arb_hypgeom_shi(arb_t res, const arb_t z, slong prec)
{

4
doc/source/arb_hypgeom.rst
View file

@ -422,7 +422,9 @@ Exponential and trigonometric integrals
Computes the sine integral of the power series *z*,
truncated to length *len*.
.. function:: void arb_hypgeom_ci(arb_t res, const arb_t z, slong prec)
.. function:: void _arb_hypgeom_ci_asymp(arb_t res, const arb_t z, slong N, slong prec)
void _arb_hypgeom_ci_2f3(arb_t res, const arb_t z, slong N, slong wp, slong prec)
void arb_hypgeom_ci(arb_t res, const arb_t z, slong prec)
Computes the cosine integral `\operatorname{Ci}(z)`.
The result is indeterminate if `z < 0` since the value of the function would be complex.