Commit c56d7ba3 authored by Julien Puydt's avatar Julien Puydt

New upstream version 2.10.0

parent cd5e3450
......@@ -13,8 +13,8 @@ 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_modular dirichlet acb_dirichlet arb_hypgeom bernoulli hypgeom \
fmpz_extras bool_mat partitions dlog \
acb_calc acb_hypgeom acb_elliptic acb_modular dirichlet acb_dirichlet \
arb_hypgeom bernoulli hypgeom fmpz_extras bool_mat partitions dlog \
$(EXTRA_BUILD_DIRS)
TEMPLATE_DIRS =
......
......@@ -50,26 +50,8 @@ acb_init(acb_t x)
void acb_clear(acb_t x);
ACB_INLINE acb_ptr
_acb_vec_init(slong n)
{
slong i;
acb_ptr v = (acb_ptr) flint_malloc(sizeof(acb_struct) * n);
for (i = 0; i < n; i++)
acb_init(v + i);
return v;
}
ACB_INLINE void
_acb_vec_clear(acb_ptr v, slong n)
{
slong i;
for (i = 0; i < n; i++)
acb_clear(v + i);
flint_free(v);
}
acb_ptr _acb_vec_init(slong n);
void _acb_vec_clear(acb_ptr v, slong n);
ACB_INLINE arb_ptr acb_real_ptr(acb_t z) { return acb_realref(z); }
ACB_INLINE arb_ptr acb_imag_ptr(acb_t z) { return acb_imagref(z); }
......@@ -110,6 +92,12 @@ acb_is_int(const acb_t z)
return arb_is_zero(acb_imagref(z)) && arb_is_int(acb_realref(z));
}
ACB_INLINE int
acb_is_int_2exp_si(const acb_t z, slong e)
{
return arb_is_zero(acb_imagref(z)) && arb_is_int_2exp_si(acb_realref(z), e);
}
ACB_INLINE void
acb_zero(acb_t z)
{
......
......@@ -79,7 +79,7 @@ acb_div(acb_t z, const acb_t x, const acb_t y, slong prec)
}
else
{
if (prec > 256 && acb_bits(y) <= prec / 2)
if (prec > 256 && acb_bits(y) <= prec / 2 && acb_is_exact(y))
{
arb_t t, u, v;
......
/*
Copyright (C) 2013 Fredrik Johansson
Copyright (C) 2017 Fredrik Johansson
This file is part of Arb.
......@@ -11,44 +11,189 @@
#include "acb.h"
static void
_arb_arf_div_rounded_den(arb_t res, const arf_t x, const arf_t y, int yinexact, slong prec)
{
int inexact = arf_div(arb_midref(res), x, y, prec, ARB_RND);
if (yinexact && !arf_is_special(arb_midref(res)))
arf_mag_set_ulp(arb_radref(res), arb_midref(res), prec - 1);
else if (inexact)
arf_mag_set_ulp(arb_radref(res), arb_midref(res), prec);
else
mag_zero(arb_radref(res));
}
static void
_arb_arf_div_rounded_den_add_err(arb_t res, const arf_t x, const arf_t y, int yinexact, slong prec)
{
int inexact = arf_div(arb_midref(res), x, y, prec, ARB_RND);
if (yinexact && !arf_is_special(arb_midref(res)))
arf_mag_add_ulp(arb_radref(res), arb_radref(res), arb_midref(res), prec - 1);
else if (inexact)
arf_mag_add_ulp(arb_radref(res), arb_radref(res), arb_midref(res), prec);
}
void
acb_inv(acb_t z, const acb_t x, slong prec)
acb_inv(acb_t res, const acb_t z, slong prec)
{
#define a acb_realref(x)
#define b acb_imagref(x)
#define c acb_realref(z)
#define d acb_imagref(z)
mag_t am, bm;
slong hprec;
#define a arb_midref(acb_realref(z))
#define b arb_midref(acb_imagref(z))
#define x arb_radref(acb_realref(z))
#define y arb_radref(acb_imagref(z))
/* choose precision for the floating-point approximation of a^2+b^2 so
that the double rounding result in less than
2 ulp error; also use at least MAG_BITS bits since the
value will be recycled for error bounds */
hprec = FLINT_MAX(prec + 3, MAG_BITS);
if (arb_is_zero(acb_imagref(z)))
{
arb_inv(acb_realref(res), acb_realref(z), prec);
arb_zero(acb_imagref(res));
return;
}
if (arb_is_zero(acb_realref(z)))
{
arb_inv(acb_imagref(res), acb_imagref(z), prec);
arb_neg(acb_imagref(res), acb_imagref(res));
arb_zero(acb_realref(res));
return;
}
if (!acb_is_finite(z))
{
acb_indeterminate(res);
return;
}
if (arb_is_zero(b))
if (mag_is_zero(x) && mag_is_zero(y))
{
arb_inv(c, a, prec);
arb_zero(d);
int inexact;
arf_t a2b2;
arf_init(a2b2);
inexact = arf_sosq(a2b2, a, b, hprec, ARF_RND_DOWN);
if (arf_is_special(a2b2))
{
acb_indeterminate(res);
}
else
{
_arb_arf_div_rounded_den(acb_realref(res), a, a2b2, inexact, prec);
_arb_arf_div_rounded_den(acb_imagref(res), b, a2b2, inexact, prec);
arf_neg(arb_midref(acb_imagref(res)), arb_midref(acb_imagref(res)));
}
arf_clear(a2b2);
return;
}
else if (arb_is_zero(a))
mag_init(am);
mag_init(bm);
/* first bound |a|-x, |b|-y */
arb_get_mag_lower(am, acb_realref(z));
arb_get_mag_lower(bm, acb_imagref(z));
if ((mag_is_zero(am) && mag_is_zero(bm)))
{
arb_inv(d, b, prec);
arb_neg(d, d);
arb_zero(c);
acb_indeterminate(res);
}
else
{
arb_t t;
arb_init(t);
/*
The propagated error in the real part is given exactly by
(a+x')/((a+x')^2+(b+y'))^2 - a/(a^2+b^2) = P / Q,
P = [(b^2-a^2) x' - a (x'^2+y'^2 + 2y'b)]
Q = [(a^2+b^2)((a+x')^2+(b+y')^2)]
where |x'| <= x and |y'| <= y, and analogously for the imaginary part.
*/
mag_t t, u, v, w;
arf_t a2b2;
int inexact;
mag_init(t);
mag_init(u);
mag_init(v);
mag_init(w);
arf_init(a2b2);
inexact = arf_sosq(a2b2, a, b, hprec, ARF_RND_DOWN);
/* compute denominator */
/* t = (|a|-x)^2 + (|b|-x)^2 (lower bound) */
mag_mul_lower(t, am, am);
mag_mul_lower(u, bm, bm);
mag_add_lower(t, t, u);
/* u = a^2 + b^2 (lower bound) */
arf_get_mag_lower(u, a2b2);
/* t = ((|a|-x)^2 + (|b|-x)^2)(a^2 + b^2) (lower bound) */
mag_mul_lower(t, t, u);
/* compute numerator */
/* real: |a^2-b^2| x + |a| ((x^2 + y^2) + 2 |b| y)) */
/* imag: |a^2-b^2| y + |b| ((x^2 + y^2) + 2 |a| x)) */
/* am, bm = upper bounds for a, b */
arf_get_mag(am, a);
arf_get_mag(bm, b);
/* v = x^2 + y^2 */
mag_mul(v, x, x);
mag_addmul(v, y, y);
/* u = |a| ((x^2 + y^2) + 2 |b| y) */
mag_mul_2exp_si(u, bm, 1);
mag_mul(u, u, y);
mag_add(u, u, v);
mag_mul(u, u, am);
/* v = |b| ((x^2 + y^2) + 2 |a| x) */
mag_mul_2exp_si(w, am, 1);
mag_addmul(v, w, x);
mag_mul(v, v, bm);
/* w = |b^2 - a^2| (upper bound) */
if (arf_cmpabs(a, b) >= 0)
mag_mul(w, am, am);
else
mag_mul(w, bm, bm);
mag_addmul(u, w, x);
mag_addmul(v, w, y);
arb_mul(t, a, a, prec);
arb_addmul(t, b, b, prec);
mag_div(arb_radref(acb_realref(res)), u, t);
mag_div(arb_radref(acb_imagref(res)), v, t);
arb_div(c, a, t, prec);
arb_div(d, b, t, prec);
_arb_arf_div_rounded_den_add_err(acb_realref(res), a, a2b2, inexact, prec);
_arb_arf_div_rounded_den_add_err(acb_imagref(res), b, a2b2, inexact, prec);
arf_neg(arb_midref(acb_imagref(res)), arb_midref(acb_imagref(res)));
arb_neg(d, d);
mag_clear(t);
mag_clear(u);
mag_clear(v);
mag_clear(w);
arb_clear(t);
arf_clear(a2b2);
}
mag_clear(am);
mag_clear(bm);
#undef a
#undef b
#undef c
#undef d
#undef x
#undef y
}
......@@ -30,7 +30,7 @@ acb_polygamma(acb_t res, const acb_t s, const acb_t z, slong prec)
acb_gamma(u, t, prec);
acb_hurwitz_zeta(t, t, z, prec);
if (arf_is_int_2exp_si(arb_midref(acb_realref(s)), 1))
if (acb_is_int_2exp_si(s, 1))
acb_neg(t, t);
acb_mul(res, t, u, prec);
......
......@@ -16,6 +16,7 @@ acb_sqrt(acb_t y, const acb_t x, slong prec)
{
arb_t r, t, u;
slong wp;
int done;
#define a acb_realref(x)
#define b acb_imagref(x)
......@@ -64,32 +65,70 @@ acb_sqrt(acb_t y, const acb_t x, slong prec)
arb_init(t);
arb_init(u);
/* r = |a+bi| */
acb_abs(r, x, wp);
arb_add(t, r, a, wp);
if (arb_rel_accuracy_bits(t) > 8)
done = 0;
if (arf_sgn(arb_midref(a)) >= 0)
{
/* sqrt(a+bi) = sqrt((r+a)/2) + b/sqrt(2*(r+a))*i, r = |a+bi| */
arb_add(t, r, a, wp);
arb_mul_2exp_si(u, t, 1);
arb_sqrt(u, u, wp);
arb_div(d, b, u, prec);
if (arb_rel_accuracy_bits(t) > 8)
{
/* sqrt(a+bi) = sqrt((r+a)/2) + b/sqrt(2*(r+a))*i */
arb_mul_2exp_si(u, t, 1);
arb_sqrt(u, u, wp);
arb_div(d, b, u, prec);
arb_set_round(c, u, prec);
arb_mul_2exp_si(c, c, -1);
done = 1;
}
else
{
arb_sub(u, r, a, wp);
}
}
else if (!arb_contains_zero(b))
{
arb_sub(u, r, a, wp);
arb_set_round(c, u, prec);
arb_mul_2exp_si(c, c, -1);
if (arb_rel_accuracy_bits(u) > 8)
{
/* sqrt(a+bi) = |b|/sqrt(2*(r-a)) + sgn(b)*sqrt((r-a)/2)*i */
int sgn = arf_sgn(arb_midref(b));
arb_mul_2exp_si(t, u, 1);
arb_sqrt(t, t, wp);
arb_div(c, b, t, prec);
arb_abs(c, c);
arb_set_round(d, t, prec);
arb_mul_2exp_si(d, d, -1);
if (sgn < 0)
arb_neg(d, d);
done = 1;
}
else
{
arb_add(t, r, a, wp);
}
}
else
{
/*
sqrt(a+bi) = sqrt((r+a)/2) + (b/|b|)*sqrt((r-a)/2)*i
(sign)
*/
arb_add(t, r, a, wp);
arb_sub(u, r, a, wp);
}
arb_mul_2exp_si(t, t, -1);
/* t = r+a, u = r-a */
arb_sub(u, r, a, wp);
if (!done)
{
/* sqrt(a+bi) = sqrt((r+a)/2) + (b/|b|)*sqrt((r-a)/2)*i
(sign) */
arb_mul_2exp_si(t, t, -1);
arb_mul_2exp_si(u, u, -1);
arb_sqrtpos(c, t, prec);
if (arb_is_nonnegative(b))
......
/*
Copyright (C) 2012 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"
static void
acb_inv_naive(acb_t z, const acb_t x, slong prec)
{
#define a acb_realref(x)
#define b acb_imagref(x)
#define c acb_realref(z)
#define d acb_imagref(z)
if (arb_is_zero(b))
{
arb_inv(c, a, prec);
arb_zero(d);
}
else if (arb_is_zero(a))
{
arb_inv(d, b, prec);
arb_neg(d, d);
arb_zero(c);
}
else
{
arb_t t;
arb_init(t);
arb_mul(t, a, a, prec);
arb_addmul(t, b, b, prec);
arb_div(c, a, t, prec);
arb_div(d, b, t, prec);
arb_neg(d, d);
arb_clear(t);
}
#undef a
#undef b
#undef c
#undef d
}
int main()
{
slong iter;
flint_rand_t state;
flint_printf("inv....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++)
{
acb_t a, b, c, d, e, f;
arf_t t;
slong prec;
acb_init(a);
acb_init(b);
acb_init(c);
acb_init(d);
acb_init(e);
acb_init(f);
arf_init(t);
prec = 2 + n_randint(state, 1000);
acb_randtest_special(a, state, 1 + n_randint(state, 1000), 100);
acb_randtest_special(b, state, 1 + n_randint(state, 1000), 100);
acb_inv(b, a, prec);
acb_inv_naive(c, a, prec);
if (!acb_overlaps(b, c))
{
flint_printf("FAIL: overlap\n\n");
flint_printf("a = "); acb_printd(a, 30); flint_printf("\n\n");
flint_printf("b = "); acb_printd(b, 30); flint_printf("\n\n");
flint_printf("c = "); acb_printd(c, 30); flint_printf("\n\n");
abort();
}
acb_set(c, a);
acb_inv(c, c, prec);
if (!acb_equal(b, c))
{
flint_printf("FAIL: aliasing\n\n");
flint_printf("a = "); acb_printd(a, 30); flint_printf("\n\n");
flint_printf("b = "); acb_printd(b, 30); flint_printf("\n\n");
flint_printf("c = "); acb_printd(c, 30); flint_printf("\n\n");
abort();
}
acb_randtest(a, state, 1 + n_randint(state, 1000), 10);
acb_randtest(b, state, 1 + n_randint(state, 1000), 10);
acb_zero(d);
arf_set_mag(t, arb_radref(acb_realref(a)));
if (n_randint(state, 2))
arf_neg(t, t);
arf_add(arb_midref(acb_realref(d)),
arb_midref(acb_realref(a)), t, ARF_PREC_EXACT, ARF_RND_DOWN);
arf_set_mag(t, arb_radref(acb_imagref(a)));
if (n_randint(state, 2))
arf_neg(t, t);
arf_add(arb_midref(acb_imagref(d)),
arb_midref(acb_imagref(a)), t, ARF_PREC_EXACT, ARF_RND_DOWN);
acb_inv(b, a, 2 + n_randint(state, 1000));
acb_inv(d, d, 2 + n_randint(state, 1000));
if (!acb_overlaps(b, d))
{
flint_printf("FAIL: corner test\n\n");
flint_printf("a = "); acb_printd(a, 30); flint_printf("\n\n");
flint_printf("b = "); acb_printd(b, 30); flint_printf("\n\n");
flint_printf("d = "); acb_printd(d, 30); flint_printf("\n\n");
abort();
}
acb_clear(a);
acb_clear(b);
acb_clear(c);
acb_clear(d);
acb_clear(e);
acb_clear(f);
arf_clear(t);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}
/*
Copyright (C) 2014 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_vec_clear(acb_ptr v, slong n)
{
slong i;
for (i = 0; i < n; i++)
acb_clear(v + i);
flint_free(v);
}
/*
Copyright (C) 2014 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"
acb_ptr
_acb_vec_init(slong n)
{
slong i;
acb_ptr v = (acb_ptr) flint_malloc(sizeof(acb_struct) * n);
for (i = 0; i < n; i++)
acb_init(v + i);
return v;
}
/*
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/>.
*/
#ifndef ACB_ELLIPTIC_H
#define ACB_ELLIPTIC_H
#include <stdio.h>
#include "acb.h"
#include "acb_poly.h"
#ifdef __cplusplus
extern "C" {
#endif
void acb_elliptic_k(acb_t k, const acb_t m, slong prec);
void acb_elliptic_k_jet(acb_ptr w, const acb_t m, slong len, slong prec);
void _acb_elliptic_k_series(acb_ptr res, acb_srcptr m, slong zlen, slong len, slong prec);
void acb_elliptic_k_series(acb_poly_t res, const acb_poly_t m, slong len, slong prec);
void acb_elliptic_e(acb_t res, const acb_t m, slong prec);
void acb_elliptic_rf(acb_t res, const acb_t x, const acb_t y, const acb_t z, int flags, slong prec);
void acb_elliptic_rj(acb_t res, const acb_t x, const acb_t y, const acb_t z, const acb_t p, int flags, slong prec);
void acb_elliptic_rg(acb_t res, const acb_t x, const acb_t y, const acb_t z, int flags, slong prec);
void acb_elliptic_rc1(acb_t res, const acb_t x, slong prec);
void acb_elliptic_f(acb_t res, const acb_t phi, const acb_t m, int times_pi, slong prec);
void acb_elliptic_e_inc(acb_t res, const acb_t phi, const acb_t m, int times_pi, slong prec);
void acb_elliptic_pi(acb_t r, const acb_t n, const acb_t m, slong prec);
void acb_elliptic_pi_inc(acb_t res, const acb_t n, const acb_t phi, const acb_t m, int times_pi, slong prec);
void acb_elliptic_p(acb_t r, const acb_t z, const acb_t tau, slong prec);
void acb_elliptic_p_jet(acb_ptr r, const acb_t z, const acb_t tau, slong len, slong prec);
void _acb_elliptic_p_series(acb_ptr res, acb_srcptr z, slong zlen, const acb_t tau, slong len, slong prec);
void acb_elliptic_p_series(acb_poly_t res, const acb_poly_t z, const acb_t tau, slong len, slong prec);
void acb_elliptic_zeta(acb_t res, const acb_t z, const acb_t tau, slong prec);
void acb_elliptic_sigma(acb_t res, const acb_t z, const acb_t tau, slong prec);
void acb_elliptic_roots(acb_t e1, acb_t e2, acb_t e3, const acb_t tau, slong prec);
void acb_elliptic_invariants(acb_t g2, acb_t g3, const acb_t tau, slong prec);
void acb_elliptic_inv_p(acb_t res, const acb_t z, const acb_t tau, slong prec);
#ifdef __cplusplus
}
#endif
#endif
/*
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_elliptic.h"
void
acb_elliptic_e(acb_t res, const acb_t m, slong prec)
{
if (acb_is_zero(m))
{
acb_const_pi(res, prec);
acb_mul_2exp_si(res, res, -1);
}
else if (acb_is_one(m))
{
acb_one(res);
}
else
{
acb_struct t[2];
acb_init(t + 0);
acb_init(t + 1);
acb_elliptic_k_jet(t, m, 2, prec);
acb_mul(t + 1, t + 1, m, prec);
acb_mul_2exp_si(t + 1, t + 1, 1);
acb_add(t, t, t + 1, prec);
acb_sub_ui(t + 1, m, 1, prec);
acb_mul(res, t, t + 1, prec);
acb_neg(res, res);
acb_clear(t + 0