Commit 8c3c7c59 authored by Julien Puydt's avatar Julien Puydt

New upstream version 2.12.0

parent e125cfc2
......@@ -17,6 +17,7 @@ compiler:
script:
- if [[ "${TRAVIS_OS_NAME}" == "osx" ]] && [[ "${CC}" == "gcc" ]]; then
brew update;
brew install gcc48;
export CC=gcc-4.8;
export CXX=g++-4.8;
......
......@@ -13,7 +13,7 @@ QUIET_AR = @echo ' ' AR ' ' $@;
AT=@
BUILD_DIRS = fmpr arf mag arb arb_mat arb_poly arb_calc acb acb_mat acb_poly \
acb_calc acb_hypgeom acb_elliptic acb_modular dirichlet acb_dirichlet \
acb_dft acb_calc acb_hypgeom acb_elliptic acb_modular dirichlet acb_dirichlet \
arb_hypgeom bernoulli hypgeom fmpz_extras bool_mat partitions dlog \
arb_fmpz_poly \
$(EXTRA_BUILD_DIRS)
......
......@@ -86,9 +86,8 @@ including methods for composition, reversion, product trees,
multipoint evaluation and interpolation, complex root isolation,
and transcendental functions of power series.
Arb has partial support for automatic differentiation (AD), and includes
rudimentary functionality for rigorous calculus based on AD
(including real root isolation and complex integration).
Other features include root isolation for real functions, rigorous numerical
integration of complex functions, and discrete Fourier transforms (DFTs).
## Special functions
......@@ -98,10 +97,10 @@ Riemann zeta and Hurwitz zeta function, Dirichlet L-functions, polylogarithm,
error function, Gauss hypergeometric function 2F1, confluent
hypergeometric functions, Bessel functions, Airy functions,
Legendre functions and other orthogonal polynomials,
exponential and trigonometric integrals, incomplete gamma function,
Jacobi theta functions, modular functions, Weierstrass elliptic function,
complete elliptic integrals, arithmetic-geometric mean,
Bernoulli numbers, partition function, Barnes G-function.
exponential and trigonometric integrals, incomplete gamma and beta functions,
Jacobi theta functions, modular functions, Weierstrass elliptic functions,
complete and incomplete elliptic integrals, arithmetic-geometric mean,
Bernoulli numbers, partition function, Barnes G-function, Lambert W function.
## Speed
......
......@@ -734,6 +734,24 @@ acb_coth(acb_t y, const acb_t x, slong prec)
acb_mul_onei(y, y);
}
void acb_sech(acb_t r, const acb_t z, slong prec);
void acb_csch(acb_t r, const acb_t z, slong prec);
ACB_INLINE void
acb_sec(acb_t y, const acb_t x, slong prec)
{
acb_mul_onei(y, x);
acb_sech(y, y, prec);
}
ACB_INLINE void
acb_csc(acb_t y, const acb_t x, slong prec)
{
acb_mul_onei(y, x);
acb_csch(y, y, prec);
acb_mul_onei(y, y);
}
void acb_sin_pi(acb_t r, const acb_t z, slong prec);
void acb_cos_pi(acb_t r, const acb_t z, slong prec);
void acb_sin_cos_pi(acb_t s, acb_t c, const acb_t z, slong prec);
......@@ -1142,7 +1160,7 @@ void _acb_vec_sort_pretty(acb_ptr vec, slong len);
/* roots of unity */
void acb_unit_root(acb_t res, ulong order, slong prec);
void _acb_vec_unit_roots(acb_ptr z, slong len, slong prec);
void _acb_vec_unit_roots(acb_ptr z, slong order, slong len, slong prec);
ACB_INLINE slong
acb_allocated_bytes(const acb_t x)
......
......@@ -73,13 +73,48 @@ sqrtmul(acb_t c, const acb_t a, const acb_t b, slong prec)
}
}
static void
acb_agm_close_taylor(acb_t res, acb_t z, acb_t z2,
const acb_t aplusb, const acb_t aminusb,
const mag_t err, slong prec)
{
acb_div(z, aminusb, aplusb, prec);
acb_sqr(z, z, prec);
acb_sqr(z2, z, prec);
acb_mul_si(res, z2, -469, prec);
acb_addmul_si(res, z, -704, prec);
acb_mul(res, res, z2, prec);
acb_addmul_si(res, z2, -1280, prec);
acb_mul_2exp_si(z, z, 12);
acb_sub(res, res, z, prec);
acb_add_ui(res, res, 16384, prec);
acb_mul_2exp_si(res, res, -15);
acb_add_error_mag(res, err);
acb_mul(res, res, aplusb, prec);
}
void
acb_agm1_basecase(acb_t res, const acb_t z, slong prec)
{
acb_t a, b, t, u;
mag_t err;
mag_t err, err2;
int isreal;
isreal = acb_is_real(z) && arb_is_nonnegative(acb_realref(z));
if (isreal)
{
acb_init(a);
acb_one(a);
arb_agm(acb_realref(res), acb_realref(a), acb_realref(z), prec);
arb_zero(acb_imagref(res));
acb_clear(a);
return;
}
if (acb_is_zero(z))
{
acb_zero(res);
......@@ -107,20 +142,43 @@ acb_agm1_basecase(acb_t res, const acb_t z, slong prec)
return;
}
isreal = acb_is_real(z) && arb_is_nonnegative(acb_realref(z));
acb_init(a);
acb_init(b);
acb_init(t);
acb_init(u);
mag_init(err);
mag_init(err2);
acb_one(a);
acb_set_round(b, z, prec);
while (!acb_overlaps(a, b))
while (1)
{
acb_sub(u, a, b, prec);
if (acb_contains_zero(u))
{
/* Dupont's thesis, p. 87: |M(z) - a_n| <= |a_n - b_n| */
acb_set(res, a);
acb_get_mag(err, u);
acb_add_error_mag(res, err);
break;
}
acb_add(t, a, b, prec);
acb_get_mag(err, u);
acb_get_mag_lower(err2, t);
mag_div(err, err, err2);
mag_geom_series(err, err, 10);
mag_mul_2exp_si(err, err, -6);
if (mag_cmp_2exp_si(err, -prec) < 0)
{
acb_agm_close_taylor(res, a, b, t, u, err, prec);
break;
}
acb_mul_2exp_si(t, t, -1);
sqrtmul(u, a, b, prec);
......@@ -129,23 +187,12 @@ acb_agm1_basecase(acb_t res, const acb_t z, slong prec)
acb_swap(u, b);
}
/* Dupont's thesis, p. 87:
|M(z) - a_n| <= |a_n - b_n| */
acb_sub(t, a, b, prec);
acb_get_mag(err, t);
if (isreal)
arb_add_error_mag(acb_realref(a), err);
else
acb_add_error_mag(a, err);
acb_set(res, a);
acb_clear(a);
acb_clear(b);
acb_clear(t);
acb_clear(u);
mag_clear(err);
mag_clear(err2);
}
/*
......@@ -155,9 +202,10 @@ acb_agm1_basecase(acb_t res, const acb_t z, slong prec)
void
acb_agm1_deriv_diff(acb_t Mz, acb_t Mzp, const acb_t z, slong prec)
{
mag_t err, t;
mag_t err, t, C;
fmpz_t rexp, hexp;
slong wp;
acb_t u, v;
slong wp, qexp;
int isreal;
if (!acb_is_exact(z) || !acb_is_finite(z) ||
......@@ -176,9 +224,13 @@ acb_agm1_deriv_diff(acb_t Mz, acb_t Mzp, const acb_t z, slong prec)
D = 1/r,
and 0 < r < |z|
M(z+h) - M(z)
|------------- - M'(z)| <= C D^2 h / (1 - D h)
h
M(z+h) - M(z-h)
|--------------- - M'(z)| <= D^3 h^2 / (1 - D h)
2h
M(z+h) + M(z-h)
|--------------- - M(z)| <= D^2 h^2 / (1 - D h)
2
h D < 1.
*/
......@@ -187,54 +239,88 @@ acb_agm1_deriv_diff(acb_t Mz, acb_t Mzp, const acb_t z, slong prec)
fmpz_init(rexp);
mag_init(err);
mag_init(t);
mag_init(C);
acb_init(u);
acb_init(v);
/* choose r = 2^rexp such that r < |z| */
acb_get_mag_lower(t, z);
fmpz_sub_ui(rexp, MAG_EXPREF(t), 2);
/* Choose h = 2^hexp with hexp = rexp - (prec + 5).
D = 2^-rexp
C D^2 h / (1 - D h) <= C * 2^(-5-prec-rexp+1)
/* Choose h = r/q = 2^hexp = 2^(rexp-qexp)
with qexp = floor(prec/2) + 5
D = 1/r = 2^-rexp
f(z) error <= C D^2 h^2 / (1-Dh)
f'(z) error <= C D^3 h^2 / (1-Dh)
1/(1-Dh) < 2, hence:
f(z) error < 2 C D^2 h^2 = C 2^(1-2*qexp)
f'(z) error < 2 C D^3 h^2 = C 2^(1-rexp-2*qexp)
*/
/* err = C = max(1, |z| + r) */
acb_get_mag(err, z);
/* C = max(1, |z| + r) */
acb_get_mag(C, z);
mag_one(t);
mag_mul_2exp_fmpz(t, t, rexp);
mag_add(err, err, t);
mag_add(C, C, t);
mag_one(t);
mag_max(err, err, t);
mag_max(C, C, t);
qexp = prec / 2 + 5;
/*
if (fmpz_sgn(rexp) < 0)
qexp += fmpz_bits(rexp);
*/
/* multiply by 2^(-5-prec-rexp+1) (use hexp as temp) */
fmpz_set_si(hexp, 1 - 5 - prec);
fmpz_sub(hexp, hexp, rexp);
mag_mul_2exp_fmpz(err, err, hexp);
/* compute h = 2^hexp */
fmpz_sub_ui(hexp, rexp, qexp);
/* choose h = 2^hexp */
fmpz_sub_ui(hexp, rexp, prec + 5);
/* compute finite differences */
wp = prec + qexp + 5;
/* compute finite difference */
wp = 2 * prec + 10;
acb_agm1_basecase(Mz, z, wp);
acb_one(Mzp);
acb_mul_2exp_fmpz(Mzp, Mzp, hexp);
acb_add(Mzp, Mzp, z, wp);
acb_agm1_basecase(Mzp, Mzp, wp);
acb_sub(Mzp, Mzp, Mz, prec);
acb_one(u);
acb_mul_2exp_fmpz(u, u, hexp);
acb_add(u, z, u, wp);
acb_agm1_basecase(u, u, wp);
acb_one(v);
acb_mul_2exp_fmpz(v, v, hexp);
acb_sub(v, z, v, wp);
acb_agm1_basecase(v, v, wp);
acb_add(Mz, u, v, prec);
acb_sub(Mzp, u, v, prec);
acb_mul_2exp_si(Mz, Mz, -1);
acb_mul_2exp_si(Mzp, Mzp, -1);
fmpz_neg(hexp, hexp);
acb_mul_2exp_fmpz(Mzp, Mzp, hexp);
/* add error */
mag_mul_2exp_si(err, C, 1 - 2 * qexp);
if (isreal)
arb_add_error_mag(acb_realref(Mz), err);
else
acb_add_error_mag(Mz, err);
fmpz_neg(rexp, rexp);
mag_mul_2exp_fmpz(err, err, rexp);
if (isreal)
arb_add_error_mag(acb_realref(Mzp), err);
else
acb_add_error_mag(Mzp, err);
acb_set_round(Mz, Mz, prec);
fmpz_clear(hexp);
fmpz_clear(rexp);
mag_clear(err);
mag_clear(t);
mag_clear(C);
acb_clear(u);
acb_clear(v);
}
/*
......@@ -369,7 +455,13 @@ acb_agm1_deriv(acb_t Mz, acb_t Mzp, const acb_t z, slong prec)
if (arf_sgn(arb_midref(acb_realref(z))) >= 0)
{
acb_agm1_deriv_right(Mz, Mzp, z, prec);
if (acb_is_one(z))
{
acb_one(Mz);
acb_mul_2exp_si(Mzp, Mz, -1);
}
else
acb_agm1_deriv_right(Mz, Mzp, z, prec);
}
else
{
......
......@@ -61,18 +61,49 @@ acb_atan(acb_t r, const acb_t z, slong prec)
acb_log1p(t, t, prec);
acb_log1p(u, u, prec);
acb_sub(t, t, u, prec);
acb_mul_onei(t, t);
acb_mul_2exp_si(r, t, -1);
}
else
else if (acb_is_exact(z))
{
acb_onei(t);
acb_sub(t, t, z, prec);
acb_div(t, z, t, prec);
acb_mul_2exp_si(t, t, 1);
acb_log1p(t, t, prec);
acb_mul_onei(t, t);
acb_mul_2exp_si(r, t, -1);
}
else
{
mag_t err, err2;
mag_init(err);
mag_init(err2);
acb_mul_onei(t, t);
acb_mul_2exp_si(r, t, -1);
/* atan'(z) = 1/(1+z^2) = 1/((z+i)(z-i)) */
acb_onei(t);
acb_add(t, z, t, prec);
acb_get_mag_lower(err, t);
acb_onei(t);
acb_sub(t, z, t, prec);
acb_get_mag_lower(err2, t);
mag_mul_lower(err, err, err2);
mag_hypot(err2, arb_radref(acb_realref(z)), arb_radref(acb_imagref(z)));
mag_div(err, err2, err);
arf_set(arb_midref(acb_realref(t)), arb_midref(acb_realref(z)));
arf_set(arb_midref(acb_imagref(t)), arb_midref(acb_imagref(z)));
mag_zero(arb_radref(acb_realref(t)));
mag_zero(arb_radref(acb_imagref(t)));
acb_atan(r, t, prec);
acb_add_error_mag(r, err);
mag_clear(err);
mag_clear(err2);
}
acb_clear(t);
acb_clear(u);
......
/*
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 "acb.h"
void
acb_csch(acb_t res, const acb_t z, slong prec)
{
if (acb_contains_zero(z) || !acb_is_finite(z))
{
acb_indeterminate(res);
}
else if (arb_is_zero(acb_imagref(z)))
{
arb_csch(acb_realref(res), acb_realref(z), prec);
arb_zero(acb_imagref(res));
}
else if (arb_is_zero(acb_realref(z)))
{
arb_csc(acb_imagref(res), acb_imagref(z), prec);
arb_neg(acb_imagref(res), acb_imagref(res));
arb_zero(acb_realref(res));
}
else
{
if (arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 0) > 0)
{
acb_t t;
acb_init(t);
if (arf_sgn(arb_midref(acb_realref(z))) > 0)
{
acb_neg(t, z);
acb_exp(t, t, prec + 4);
acb_mul(res, t, t, prec + 4);
acb_sub_ui(res, res, 1, prec + 4);
acb_div(res, t, res, prec);
acb_neg(res, res);
}
else
{
acb_exp(t, z, prec + 4);
acb_mul(res, t, t, prec + 4);
acb_sub_ui(res, res, 1, prec + 4);
acb_div(res, t, res, prec);
}
acb_mul_2exp_si(res, res, 1);
acb_clear(t);
}
else
{
acb_sinh(res, z, prec + 4);
acb_inv(res, res, prec);
}
}
}
......@@ -56,8 +56,18 @@ acb_pow_fmpz_binexp(acb_t y, const acb_t b, const fmpz_t e, slong prec)
fmpz_t f;
fmpz_init(f);
fmpz_neg(f, e);
acb_pow_fmpz_binexp(y, b, f, prec + 2);
acb_inv(y, y, prec);
if (acb_is_exact(b))
{
acb_pow_fmpz_binexp(y, b, f, prec + 2);
acb_inv(y, y, prec);
}
else
{
acb_inv(y, b, prec + fmpz_bits(e) + 2);
acb_pow_fmpz_binexp(y, y, f, prec);
}
fmpz_clear(f);
return;
}
......@@ -180,8 +190,17 @@ acb_pow_arb(acb_t z, const acb_t x, const arb_t y, slong prec)
else
{
arf_get_fmpz_fixed_si(e, ymid, -1);
acb_sqrt(z, x, prec + fmpz_bits(e));
acb_pow_fmpz_binexp(z, z, e, prec);
if (fmpz_sgn(e) >= 0)
{
acb_sqrt(z, x, prec + fmpz_bits(e));
acb_pow_fmpz_binexp(z, z, e, prec);
}
else
{
fmpz_neg(e, e);
acb_rsqrt(z, x, prec + fmpz_bits(e));
acb_pow_fmpz_binexp(z, z, e, prec);
}
}
fmpz_clear(e);
......
/*
Copyright (C) 2017 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 "acb.h"
void
acb_sech(acb_t res, const acb_t z, slong prec)
{
if (arb_is_zero(acb_imagref(z)))
{
arb_sech(acb_realref(res), acb_realref(z), prec);
arb_zero(acb_imagref(res));
}
else if (arb_is_zero(acb_realref(z)))
{
arb_sec(acb_realref(res), acb_imagref(z), prec);
arb_zero(acb_imagref(res));
}
else
{
if (arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 0) > 0)
{
acb_t t;
acb_init(t);
if (arf_sgn(arb_midref(acb_realref(z))) > 0)
{
acb_neg(t, z);
acb_exp(t, t, prec + 4);
}
else
{
acb_exp(t, z, prec + 4);
}
acb_mul(res, t, t, prec + 4);
acb_add_ui(res, res, 1, prec + 4);
acb_div(res, t, res, prec);
acb_mul_2exp_si(res, res, 1);
acb_clear(t);
}
else
{
acb_cosh(res, z, prec + 4);
acb_inv(res, res, prec);
}
}
}
/*
Copyright (C) 2017 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 "acb.h"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("csc....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
{
acb_t x, a, b;
slong prec1, prec2;
prec1 = 2 + n_randint(state, 200);
prec2 = prec1 + 30;
acb_init(x);
acb_init(a);
acb_init(b);
acb_randtest_special(x, state, 1 + n_randint(state, 300), 100);
acb_randtest_special(a, state, 1 + n_randint(state, 300), 100);
acb_randtest_special(b, state, 1 + n_randint(state, 300), 100);
if (n_randint(state, 2))
{
acb_csc(a, x, prec1);
}
else
{
acb_set(a, x);
acb_csc(a, a, prec1);
}
acb_sin(b, x, prec2);
acb_inv(b, b, prec2);
/* check consistency */
if (!acb_overlaps(a, b))
{
flint_printf("FAIL: overlap\n\n");
flint_printf("x = "); acb_printn(x, 20, 0); flint_printf("\n\n");
flint_printf("a = "); acb_printn(a, 20, 0); flint_printf("\n\n");
flint_printf("b = "); acb_printn(b, 20, 0); flint_printf("\n\n");
flint_abort();
}
acb_clear(x);
acb_clear(a);
acb_clear(b);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}
/*
Copyright (C) 2017 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 "acb.h"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("csch....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
{
acb_t x, a, b;
slong prec1, prec2;
prec1 = 2 + n_randint(state, 200);
prec2 = prec1 + 30;
acb_init(x);
acb_init(a);
acb_init(b);
acb_randtest_special(x, state, 1 + n_randint(state, 300), 100);
acb_randtest_special(a, state, 1 + n_randint(state, 300), 100);