mirror of
https://github.com/vale981/arb
synced 2025-03-05 09:21:38 -05:00
implement the complex arithmetic-geometric mean (code is still messy)
This commit is contained in:
Fredrik Johansson
parent
aa14f864a3
commit
c6395d0bd6
4 changed files with 759 additions and 0 deletions
3
acb.h
3
acb.h
|
|
@ -616,6 +616,9 @@ void acb_hurwitz_zeta(acb_t z, const acb_t s, const acb_t a, long prec);
|
|||
void acb_polylog(acb_t w, const acb_t s, const acb_t z, long prec);
|
||||
void acb_polylog_si(acb_t w, long s, const acb_t z, long prec);
|
||||
|
||||
void acb_agm1(acb_t m, const acb_t z, long prec);
|
||||
void acb_agm1_cpx(acb_ptr m, const acb_t z, long len, long prec);
|
||||
|
||||
/*
|
||||
TBD
|
||||
|
||||
|
|
|
|||
537
acb/agm1.c
Normal file
537
acb/agm1.c
Normal file
|
|
@ -0,0 +1,537 @@
|
|||
/*=============================================================================
|
||||
|
||||
This file is part of ARB.
|
||||
|
||||
ARB is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
ARB is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with ARB; if not, write to the Free Software
|
||||
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
|
||||
|
||||
=============================================================================*/
|
||||
/******************************************************************************
|
||||
|
||||
Copyright (C) 2014 Fredrik Johansson
|
||||
|
||||
******************************************************************************/
|
||||
|
||||
#include "acb.h"
|
||||
#include "acb_poly.h"
|
||||
|
||||
void
|
||||
mag_hypot(mag_t z, const mag_t x, const mag_t y)
|
||||
{
|
||||
if (mag_is_zero(y))
|
||||
{
|
||||
mag_set(z, x);
|
||||
}
|
||||
else if (mag_is_zero(x))
|
||||
{
|
||||
mag_set(z, y);
|
||||
}
|
||||
else
|
||||
{
|
||||
mag_t t;
|
||||
mag_init(t);
|
||||
mag_mul(t, x, x);
|
||||
mag_addmul(t, y, y);
|
||||
mag_sqrt(z, t);
|
||||
mag_clear(t);
|
||||
}
|
||||
}
|
||||
|
||||
static __inline__ void
|
||||
acb_mul_2exp_fmpz(acb_t z, const acb_t x, const fmpz_t c)
|
||||
{
|
||||
arb_mul_2exp_fmpz(acb_realref(z), acb_realref(x), c);
|
||||
arb_mul_2exp_fmpz(acb_imagref(z), acb_imagref(x), c);
|
||||
}
|
||||
|
||||
/* Checks that |arg(z)| <= 3 pi / 4 */
|
||||
static int
|
||||
acb_check_arg(const acb_t z)
|
||||
{
|
||||
mag_t re, im;
|
||||
int res;
|
||||
|
||||
if (!arb_contains_negative(acb_realref(z)))
|
||||
return 1;
|
||||
|
||||
mag_init(re);
|
||||
mag_init(im);
|
||||
|
||||
arb_get_mag(re, acb_realref(z));
|
||||
arb_get_mag_lower(im, acb_imagref(z));
|
||||
|
||||
res = mag_cmp(re, im) < 0;
|
||||
|
||||
mag_clear(re);
|
||||
mag_clear(im);
|
||||
|
||||
return res;
|
||||
}
|
||||
|
||||
static void
|
||||
sqrtmul(acb_t c, const acb_t a, const acb_t b, long prec)
|
||||
{
|
||||
if (arb_is_positive(acb_realref(a)) &&
|
||||
arb_is_positive(acb_realref(b)))
|
||||
{
|
||||
acb_mul(c, a, b, prec);
|
||||
acb_sqrt(c, c, prec);
|
||||
}
|
||||
else if (arb_is_nonnegative(acb_imagref(a)) &&
|
||||
arb_is_nonnegative(acb_imagref(b)))
|
||||
{
|
||||
acb_mul(c, a, b, prec);
|
||||
acb_neg(c, c);
|
||||
acb_sqrt(c, c, prec);
|
||||
acb_mul_onei(c, c);
|
||||
}
|
||||
else if (arb_is_nonpositive(acb_imagref(a)) &&
|
||||
arb_is_nonpositive(acb_imagref(b)))
|
||||
{
|
||||
acb_mul(c, a, b, prec);
|
||||
acb_neg(c, c);
|
||||
acb_sqrt(c, c, prec);
|
||||
acb_mul_onei(c, c);
|
||||
acb_neg(c, c);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_t d;
|
||||
acb_init(d);
|
||||
acb_sqrt(c, a, prec);
|
||||
acb_sqrt(d, b, prec);
|
||||
acb_mul(c, c, d, prec);
|
||||
acb_clear(d);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
acb_agm1_basecase(acb_t res, const acb_t z, long prec)
|
||||
{
|
||||
acb_t a, b, t, u;
|
||||
mag_t err;
|
||||
int isreal;
|
||||
|
||||
if (acb_is_zero(z))
|
||||
{
|
||||
acb_zero(res);
|
||||
return;
|
||||
}
|
||||
|
||||
if (acb_is_one(z))
|
||||
{
|
||||
acb_one(res);
|
||||
return;
|
||||
}
|
||||
|
||||
if (!acb_check_arg(z))
|
||||
{
|
||||
mag_t one;
|
||||
mag_init(one);
|
||||
mag_init(err);
|
||||
mag_one(one);
|
||||
acb_get_mag(err, z);
|
||||
mag_max(err, err, one);
|
||||
acb_zero(res);
|
||||
acb_add_error_mag(res, err);
|
||||
mag_clear(err);
|
||||
mag_clear(one);
|
||||
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);
|
||||
|
||||
acb_one(a);
|
||||
acb_set_round(b, z, prec);
|
||||
|
||||
while (!acb_overlaps(a, b))
|
||||
{
|
||||
acb_add(t, a, b, prec);
|
||||
acb_mul_2exp_si(t, t, -1);
|
||||
|
||||
sqrtmul(u, a, b, prec);
|
||||
|
||||
acb_swap(t, a);
|
||||
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);
|
||||
}
|
||||
|
||||
/*
|
||||
Computes (M(z), M'(z)) using a finite difference.
|
||||
Assumes z exact, |arg(z)| <= 3 pi / 4.
|
||||
*/
|
||||
void
|
||||
acb_agm1_deriv_diff(acb_t Mz, acb_t Mzp, const acb_t z, long prec)
|
||||
{
|
||||
mag_t err, t;
|
||||
fmpz_t rexp, hexp;
|
||||
long wp;
|
||||
int isreal;
|
||||
|
||||
if (!acb_is_exact(z) || !acb_is_finite(z) ||
|
||||
acb_is_zero(z) || !acb_check_arg(z))
|
||||
{
|
||||
acb_indeterminate(Mz);
|
||||
acb_indeterminate(Mzp);
|
||||
return;
|
||||
}
|
||||
|
||||
isreal = acb_is_real(z) && arb_is_nonnegative(acb_realref(z));
|
||||
|
||||
/*
|
||||
|M^(k)(z) / k!| <= C * D^k where
|
||||
C = max(1, |z| + r),
|
||||
D = 1/r,
|
||||
and 0 < r < |z|
|
||||
|
||||
M(z+h) - M(z)
|
||||
|------------- - M'(z)| <= C D^2 h / (1 - D h)
|
||||
h
|
||||
|
||||
h D < 1.
|
||||
*/
|
||||
|
||||
fmpz_init(hexp);
|
||||
fmpz_init(rexp);
|
||||
mag_init(err);
|
||||
mag_init(t);
|
||||
|
||||
/* 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)
|
||||
*/
|
||||
|
||||
/* err = C = max(1, |z| + r) */
|
||||
acb_get_mag(err, z);
|
||||
mag_one(t);
|
||||
mag_mul_2exp_fmpz(t, t, rexp);
|
||||
mag_add(err, err, t);
|
||||
mag_one(t);
|
||||
mag_max(err, err, t);
|
||||
|
||||
/* 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);
|
||||
|
||||
/* choose h = 2^hexp */
|
||||
fmpz_sub_ui(hexp, rexp, prec + 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);
|
||||
fmpz_neg(hexp, hexp);
|
||||
acb_mul_2exp_fmpz(Mzp, Mzp, hexp);
|
||||
|
||||
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);
|
||||
}
|
||||
|
||||
/*
|
||||
For input z + eps
|
||||
|
||||
First derivative bound: max(1, |z|+|eps|+r) / r
|
||||
Second derivative bound: 2 max(1, |z|+|eps|+r) / r^2
|
||||
|
||||
This is assuming that the circle at z with radius |eps| + r
|
||||
does not cross the negative half axis, which we check.
|
||||
*/
|
||||
|
||||
void
|
||||
acb_agm1_deriv_right(acb_t Mz, acb_t Mzp, const acb_t z, long prec)
|
||||
{
|
||||
if (acb_is_exact(z))
|
||||
{
|
||||
acb_agm1_deriv_diff(Mz, Mzp, z, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
if (!acb_is_finite(z) || !acb_check_arg(z))
|
||||
{
|
||||
acb_indeterminate(Mz);
|
||||
acb_indeterminate(Mzp);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_t t;
|
||||
mag_t r, eps, err, one;
|
||||
int isreal;
|
||||
|
||||
acb_init(t);
|
||||
mag_init(r);
|
||||
mag_init(err);
|
||||
mag_init(one);
|
||||
mag_init(eps);
|
||||
|
||||
isreal = acb_is_real(z) && arb_is_nonnegative(acb_realref(z));
|
||||
|
||||
mag_hypot(eps, arb_radref(acb_realref(z)), arb_radref(acb_imagref(z)));
|
||||
|
||||
/* choose r avoiding overlap with negative half axis */
|
||||
if (arf_sgn(arb_midref(acb_realref(z))) < 0)
|
||||
arb_get_mag_lower(r, acb_imagref(z));
|
||||
else
|
||||
acb_get_mag_lower(r, z);
|
||||
|
||||
mag_mul_2exp_si(r, r, -1);
|
||||
|
||||
if (mag_is_zero(r))
|
||||
{
|
||||
acb_indeterminate(Mz);
|
||||
acb_indeterminate(Mzp);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_set(t, z);
|
||||
mag_zero(arb_radref(acb_realref(t)));
|
||||
mag_zero(arb_radref(acb_imagref(t)));
|
||||
|
||||
acb_get_mag(err, z);
|
||||
mag_add(err, err, r);
|
||||
mag_add(err, err, eps);
|
||||
mag_one(one);
|
||||
mag_max(err, err, one);
|
||||
mag_mul(err, err, eps);
|
||||
|
||||
acb_agm1_deriv_diff(Mz, Mzp, t, prec);
|
||||
|
||||
mag_div(err, err, r);
|
||||
|
||||
if (isreal)
|
||||
arb_add_error_mag(acb_realref(Mz), err);
|
||||
else
|
||||
acb_add_error_mag(Mz, err);
|
||||
|
||||
mag_div(err, err, r);
|
||||
mag_mul_2exp_si(err, err, 1);
|
||||
|
||||
if (isreal)
|
||||
arb_add_error_mag(acb_realref(Mzp), err);
|
||||
else
|
||||
acb_add_error_mag(Mzp, err);
|
||||
}
|
||||
|
||||
acb_clear(t);
|
||||
mag_clear(r);
|
||||
mag_clear(err);
|
||||
mag_clear(one);
|
||||
mag_clear(eps);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
acb_agm1(acb_t m, const acb_t z, long prec)
|
||||
{
|
||||
if (arf_sgn(arb_midref(acb_realref(z))) >= 0)
|
||||
{
|
||||
acb_agm1_basecase(m, z, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
/* use M(z) = (z+1)/2 * M(2 sqrt(z) / (z+1)) */
|
||||
acb_t t;
|
||||
acb_init(t);
|
||||
acb_add_ui(t, z, 1, prec);
|
||||
acb_sqrt(m, z, prec);
|
||||
acb_div(m, m, t, prec);
|
||||
acb_mul_2exp_si(m, m, 1);
|
||||
acb_agm1_basecase(m, m, prec);
|
||||
acb_mul(m, m, t, prec);
|
||||
acb_mul_2exp_si(m, m, -1);
|
||||
acb_clear(t);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
acb_agm1_deriv(acb_t Mz, acb_t Mzp, const acb_t z, long prec)
|
||||
{
|
||||
/*
|
||||
u = 2 sqrt(z) / (1+z)
|
||||
|
||||
Mz = (1+z) M(u) / 2
|
||||
Mzp = [M(u) - (z-1) M'(u) / ((1+z) sqrt(z))] / 2
|
||||
*/
|
||||
|
||||
if (arf_sgn(arb_midref(acb_realref(z))) >= 0)
|
||||
{
|
||||
acb_agm1_deriv_right(Mz, Mzp, z, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_t t, u, zp1, zm1;
|
||||
|
||||
acb_init(t);
|
||||
acb_init(u);
|
||||
acb_init(zp1);
|
||||
acb_init(zm1);
|
||||
|
||||
acb_sqrt(t, z, prec);
|
||||
acb_add_ui(zp1, z, 1, prec);
|
||||
acb_sub_ui(zm1, z, 1, prec);
|
||||
|
||||
acb_div(u, t, zp1, prec);
|
||||
acb_mul_2exp_si(u, u, 1);
|
||||
|
||||
acb_agm1_deriv_right(Mz, Mzp, u, prec);
|
||||
|
||||
acb_mul(Mzp, Mzp, zm1, prec);
|
||||
acb_mul(t, t, zp1, prec);
|
||||
acb_div(Mzp, Mzp, t, prec);
|
||||
acb_sub(Mzp, Mz, Mzp, prec);
|
||||
acb_mul_2exp_si(Mzp, Mzp, -1);
|
||||
|
||||
acb_mul(Mz, Mz, zp1, prec);
|
||||
acb_mul_2exp_si(Mz, Mz, -1);
|
||||
|
||||
acb_clear(t);
|
||||
acb_clear(u);
|
||||
acb_clear(zp1);
|
||||
acb_clear(zm1);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
acb_agm1_cpx(acb_ptr m, const acb_t z, long len, long prec)
|
||||
{
|
||||
if (len < 1)
|
||||
return;
|
||||
|
||||
if (len == 1)
|
||||
{
|
||||
acb_agm1(m, z, prec);
|
||||
return;
|
||||
}
|
||||
|
||||
if (len == 2)
|
||||
{
|
||||
acb_agm1_deriv(m, m + 1, z, prec);
|
||||
return;
|
||||
}
|
||||
|
||||
if (len >= 3)
|
||||
{
|
||||
acb_t t, u, v;
|
||||
acb_ptr w;
|
||||
long k, n;
|
||||
|
||||
acb_init(t);
|
||||
acb_init(u);
|
||||
acb_init(v);
|
||||
w = _acb_vec_init(len);
|
||||
|
||||
acb_agm1_deriv(w, w + 1, z, prec);
|
||||
|
||||
/* invert series */
|
||||
acb_inv(w, w, prec);
|
||||
acb_mul(t, w, w, prec);
|
||||
acb_mul(w + 1, w + 1, t, prec);
|
||||
acb_neg(w + 1, w + 1);
|
||||
|
||||
if (acb_is_one(z))
|
||||
{
|
||||
for (k = 2; k < len; k++)
|
||||
{
|
||||
n = k - 2;
|
||||
|
||||
acb_mul_ui(w + k, w + n + 0, (n+1)*(n+1), prec);
|
||||
acb_addmul_ui(w + k, w + n + 1, 7+3*n*(3+n), prec);
|
||||
acb_div_ui(w + k, w + k, 2*(n+2)*(n+2), prec);
|
||||
acb_neg(w + k, w + k);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
/* t = 3z^2 - 1 */
|
||||
/* u = -1 / (z^3 - z) */
|
||||
acb_mul(t, z, z, prec);
|
||||
acb_mul(u, t, z, prec);
|
||||
acb_mul_ui(t, t, 3, prec);
|
||||
acb_sub_ui(t, t, 1, prec);
|
||||
acb_sub(u, u, z, prec);
|
||||
acb_inv(u, u, prec);
|
||||
acb_neg(u, u);
|
||||
|
||||
/* use differential equation for second derivative */
|
||||
acb_mul(w + 2, z, w + 0, prec);
|
||||
acb_addmul(w + 2, t, w + 1, prec);
|
||||
acb_mul(w + 2, w + 2, u, prec);
|
||||
acb_mul_2exp_si(w + 2, w + 2, -1);
|
||||
|
||||
/* recurrence */
|
||||
for (k = 3; k < len; k++)
|
||||
{
|
||||
n = k - 3;
|
||||
acb_mul_ui(w + k, w + n + 0, (n+1)*(n+1), prec);
|
||||
acb_mul(v, w + n + 1, z, prec);
|
||||
acb_addmul_ui(w + k, v, 7+3*n*(3+n), prec);
|
||||
acb_mul(v, w + n + 2, t, prec);
|
||||
acb_addmul_ui(w + k, v, (n+2)*(n+2), prec);
|
||||
acb_mul(w + k, w + k, u, prec);
|
||||
acb_div_ui(w + k, w + k, (n+2)*(n+3), prec);
|
||||
}
|
||||
}
|
||||
|
||||
/* invert series */
|
||||
_acb_poly_inv_series(m, w, len, len, prec);
|
||||
|
||||
acb_clear(t);
|
||||
acb_clear(u);
|
||||
acb_clear(v);
|
||||
_acb_vec_clear(w, len);
|
||||
}
|
||||
}
|
||||
|
||||
201
acb/test/t-agm1.c
Normal file
201
acb/test/t-agm1.c
Normal file
|
|
@ -0,0 +1,201 @@
|
|||
/*=============================================================================
|
||||
|
||||
This file is part of ARB.
|
||||
|
||||
ARB is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
ARB is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with ARB; if not, write to the Free Software
|
||||
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
|
||||
|
||||
=============================================================================*/
|
||||
/******************************************************************************
|
||||
|
||||
Copyright (C) 2014 Fredrik Johansson
|
||||
|
||||
******************************************************************************/
|
||||
|
||||
#include "acb.h"
|
||||
|
||||
#define EPS 1e-13
|
||||
#define NUM_DERIVS 4
|
||||
#define NUM_TESTS 56
|
||||
|
||||
const double agm_testdata[NUM_TESTS][11] = {
|
||||
{1.0, 0.0, 1.0, 0.0, 0.5, 0.0, -0.0625, 0.0, 0.03125, 0.0},
|
||||
{0.25, 0.0, 0.56075714507190064253, 0.0, 0.76633907325304843764, 0.0,
|
||||
-0.58010113691169132987, 0.0, 1.2991960360521313649, 0.0},
|
||||
{1.0, 1.0, 1.0491605287327802205, 0.47815574608816122933,
|
||||
0.44643105633549979073, -0.08043578677710866283,
|
||||
-0.015455284495904882924, 0.031976374729700173479,
|
||||
-0.005073437378084728324, -0.010673958729796444985},
|
||||
{0.0, 1.0, 0.59907011736779610372, 0.59907011736779610372,
|
||||
0.43640657965245804105, -0.16266353771533806267,
|
||||
0.031271486774469549792, 0.031271486774469549792,
|
||||
-0.0084910439043492636266, 0.022780442870120286166},
|
||||
{-1.0, 1.0, 0.18841106798868002055, 0.77800407878895828015,
|
||||
0.39320630832295102335, -0.19287323123455182026,
|
||||
0.016808115488846724979, 0.020163502567351546742,
|
||||
-0.011189063384352573532, -0.0045629824424054816356},
|
||||
{-0.25, 0.0, 0.24392673474953340413, 0.27799893427725564501,
|
||||
0.079794586546549867566, -0.32574453868033984617,
|
||||
-0.27530931370867131683, -0.60287933825809104112,
|
||||
-0.54714813521966247214, -0.97291617788393560861},
|
||||
{-2.0, 0.0, -0.42296620840880168736, 0.66126618346180476447,
|
||||
0.29655367830470777795, -0.27314834694816402295,
|
||||
0.018331225229748304014, -0.043277191398267359935,
|
||||
0.0094104572902577354423, -0.021071819981331905571},
|
||||
{-0.99999994039535522461, 0.0, 0.0069953999943208591971,
|
||||
0.08334545077423930704, 12454.444757282471906, 73670.447089584396744,
|
||||
-87884330303.90316493, -553342797575.05204396, 909784741818095264.43,
|
||||
5883796037072491540.4},
|
||||
{-1.0000000596046447754, 0.0, -0.0069954003670321135601,
|
||||
0.083345455480285282001, 12454.451634795395449, -73670.529975671147878,
|
||||
87884324285.574775284, -553342761339.97964199, 909784713383817144.98,
|
||||
-5883795855289758433.3},
|
||||
{-1.0, 5.9604644775390625e-8, 2.3677293150757997928e-9,
|
||||
0.083932595219022127005, 75242.17179390509748, -0.041726661532440829443,
|
||||
-18488.306894081923308, 563725525035.34898339, -5988414396610684162.5,
|
||||
-92537341636.835656296},
|
||||
{-1.0, -5.9604644775390625e-8, 2.3677293150757997928e-9,
|
||||
-0.083932595219022127005, 75242.17179390509748, 0.041726661532440829443,
|
||||
-18488.306894081923308, -563725525035.34898339, -5988414396610684162.5,
|
||||
92537341636.835656296},
|
||||
};
|
||||
|
||||
int main()
|
||||
{
|
||||
long iter;
|
||||
flint_rand_t state;
|
||||
|
||||
printf("agm1....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
||||
/* check particular values against table */
|
||||
{
|
||||
acb_t z, t;
|
||||
acb_ptr w1;
|
||||
long i, j, prec, cnj;
|
||||
|
||||
acb_init(z);
|
||||
acb_init(t);
|
||||
w1 = _acb_vec_init(NUM_DERIVS);
|
||||
|
||||
for (prec = 32; prec <= 512; prec *= 4)
|
||||
{
|
||||
for (i = 0; i < NUM_TESTS; i++)
|
||||
{
|
||||
for (cnj = 0; cnj < 2; cnj++)
|
||||
{
|
||||
if (cnj == 1 && agm_testdata[i][0] < 0 &&
|
||||
agm_testdata[i][1] == 0)
|
||||
continue;
|
||||
|
||||
acb_zero(z);
|
||||
arf_set_d(arb_midref(acb_realref(z)), agm_testdata[i][0]);
|
||||
arf_set_d(arb_midref(acb_imagref(z)), cnj ? -agm_testdata[i][1] : agm_testdata[i][1]);
|
||||
|
||||
acb_agm1_cpx(w1, z, NUM_DERIVS, prec);
|
||||
|
||||
for (j = 0; j < NUM_DERIVS; j++)
|
||||
{
|
||||
arf_set_d(arb_midref(acb_realref(t)), agm_testdata[i][2+2*j]);
|
||||
mag_set_d(arb_radref(acb_realref(t)), fabs(agm_testdata[i][2+2*j]) * EPS);
|
||||
arf_set_d(arb_midref(acb_imagref(t)), cnj ? -agm_testdata[i][2+2*j+1] : agm_testdata[i][2+2*j+1]);
|
||||
mag_set_d(arb_radref(acb_imagref(t)), fabs(agm_testdata[i][2+2*j+1]) * EPS);
|
||||
|
||||
if (!acb_overlaps(w1 + j, t))
|
||||
{
|
||||
printf("FAIL\n\n");
|
||||
printf("j = %ld\n\n", j);
|
||||
printf("z = "); acb_printd(z, 15); printf("\n\n");
|
||||
printf("t = "); acb_printd(t, 15); printf("\n\n");
|
||||
printf("w1 = "); acb_printd(w1 + j, 15); printf("\n\n");
|
||||
abort();
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
_acb_vec_clear(w1, NUM_DERIVS);
|
||||
acb_clear(z);
|
||||
acb_clear(t);
|
||||
}
|
||||
|
||||
/* self-consistency test */
|
||||
for (iter = 0; iter < 1000; iter++)
|
||||
{
|
||||
acb_ptr m1, m2;
|
||||
acb_t z1, z2, t;
|
||||
long i, len1, len2, prec1, prec2;
|
||||
|
||||
len1 = n_randint(state, 10);
|
||||
len2 = n_randint(state, 10);
|
||||
|
||||
prec1 = 2 + n_randint(state, 2000);
|
||||
prec2 = 2 + n_randint(state, 2000);
|
||||
|
||||
m1 = _acb_vec_init(len1);
|
||||
m2 = _acb_vec_init(len2);
|
||||
|
||||
acb_init(z1);
|
||||
acb_init(z2);
|
||||
acb_init(t);
|
||||
|
||||
acb_randtest(z1, state, prec1, 1 + n_randint(state, 100));
|
||||
|
||||
if (n_randint(state, 2))
|
||||
{
|
||||
acb_set(z2, z1);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_randtest(t, state, prec2, 1 + n_randint(state, 100));
|
||||
acb_add(z2, z1, t, prec2);
|
||||
acb_sub(z2, z2, t, prec2);
|
||||
}
|
||||
|
||||
acb_agm1_cpx(m1, z1, len1, prec1);
|
||||
acb_agm1_cpx(m2, z2, len2, prec2);
|
||||
|
||||
for (i = 0; i < FLINT_MIN(len1, len2); i++)
|
||||
{
|
||||
if (!acb_overlaps(m1 + i, m2 + i))
|
||||
{
|
||||
printf("FAIL (overlap)\n\n");
|
||||
printf("iter = %ld, i = %ld, len1 = %ld, len2 = %ld, prec1 = %ld, prec2 = %ld\n\n",
|
||||
iter, i, len1, len2, prec1, prec2);
|
||||
|
||||
printf("z1 = "); acb_printd(z1, 30); printf("\n\n");
|
||||
printf("z2 = "); acb_printd(z2, 30); printf("\n\n");
|
||||
printf("m1 = "); acb_printd(m1, 30); printf("\n\n");
|
||||
printf("m2 = "); acb_printd(m2, 30); printf("\n\n");
|
||||
abort();
|
||||
}
|
||||
}
|
||||
|
||||
_acb_vec_clear(m1, len1);
|
||||
_acb_vec_clear(m2, len2);
|
||||
|
||||
acb_clear(z1);
|
||||
acb_clear(z2);
|
||||
acb_clear(t);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
||||
|
|
@ -519,3 +519,21 @@ Polylogarithms
|
|||
.. function:: void acb_polylog_si(acb_t w, long s, const acb_t z, long prec)
|
||||
|
||||
Sets *w* to the polylogarithm `\operatorname{Li}_s(z)`.
|
||||
|
||||
Arithmetic-geometric mean
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
.. function:: void acb_agm1(acb_t m, const acb_t z, long prec)
|
||||
|
||||
Sets *m* to the arithmetic-geometric mean `M(z) = \operatorname{agm}(1,z)`,
|
||||
defined such that the function is continuous in the complex plane except for
|
||||
a branch cut along the negative half axis (where it is continuous
|
||||
from above). This corresponds to always choosing an "optimal" branch for
|
||||
the square root in the arithmetic-geometric mean iteration.
|
||||
|
||||
.. function:: void acb_agm1_cpx(acb_ptr m, const acb_t z, long len, long prec)
|
||||
|
||||
Sets the coefficients in the array *m* to the power series expansion of the
|
||||
arithmetic-geometric mean at the point *z* truncated to length *len*, i.e.
|
||||
`M(z+x) \in \mathbb{C}[[x]]`.
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue