From c6395d0bd69d680826fc375428aa03ea21c40e8b Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Tue, 16 Dec 2014 16:53:24 +0100 Subject: [PATCH] implement the complex arithmetic-geometric mean (code is still messy) --- acb.h | 3 + acb/agm1.c | 537 +++++++++++++++++++++++++++++++++++++++++++++ acb/test/t-agm1.c | 201 +++++++++++++++++ doc/source/acb.rst | 18 ++ 4 files changed, 759 insertions(+) create mode 100644 acb/agm1.c create mode 100644 acb/test/t-agm1.c diff --git a/acb.h b/acb.h index 41e85bbf..1eb33958 100644 --- a/acb.h +++ b/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 diff --git a/acb/agm1.c b/acb/agm1.c new file mode 100644 index 00000000..7e6ec033 --- /dev/null +++ b/acb/agm1.c @@ -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); + } +} + diff --git a/acb/test/t-agm1.c b/acb/test/t-agm1.c new file mode 100644 index 00000000..a62231ac --- /dev/null +++ b/acb/test/t-agm1.c @@ -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; +} + diff --git a/doc/source/acb.rst b/doc/source/acb.rst index 86636b22..e0714035 100644 --- a/doc/source/acb.rst +++ b/doc/source/acb.rst @@ -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]]`. +