Skip to content

Instantly share code, notes, and snippets.

@fredrik-johansson
Created September 13, 2018 05:09
Show Gist options
  • Save fredrik-johansson/0a1edb681144b0de977984e059ef4fb9 to your computer and use it in GitHub Desktop.
Save fredrik-johansson/0a1edb681144b0de977984e059ef4fb9 to your computer and use it in GitHub Desktop.
debug stieltjes.c
/*
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 <math.h>
#include "flint/double_extras.h"
#include "acb_dirichlet.h"
#include "acb_calc.h"
/* Bound the quadratic Taylor error term. */
static void
stieltjes_bound_quadratic_term(arb_t B, const acb_t z,
const fmpz_t n1, const acb_t alpha, slong prec)
{
acb_t t, log_t;
mag_t b, c, d;
acb_init(t);
acb_init(log_t);
mag_init(b);
mag_init(c);
mag_init(d);
/* g''(z) = (n+1)(1+1/log(t)) / (t^2 log(t)) */
/* t = alpha + iz, log_t = log(t) */
acb_mul_onei(t, z);
acb_add(t, t, alpha, prec);
acb_log(log_t, t, prec);
acb_get_mag_lower(b, t);
acb_get_mag_lower(c, log_t);
/* d = 1+1/log(t) */
mag_inv(d, c);
mag_add_ui(d, d, 1);
/* divide by t^2 log(t) */
mag_div(d, d, b);
mag_div(d, d, b);
mag_div(d, d, c);
mag_mul_fmpz(d, d, n1);
/* For Taylor remainder: 1/2 rad(z)^2 */
mag_hypot(b, arb_radref(acb_realref(z)), arb_radref(acb_imagref(z)));
mag_mul(b, b, b);
mag_mul_2exp_si(b, b, -1);
mag_mul(d, d, b);
arf_set_mag(arb_midref(B), d);
mag_zero(arb_radref(B));
acb_clear(t);
acb_clear(log_t);
mag_clear(b);
mag_clear(c);
mag_clear(d);
}
static void
stieltjes_bound_large3(arb_t B, const acb_t x, const fmpz_t n1,
const acb_t alpha, slong prec)
{
acb_t m, t, log_t, u, v;
acb_t g0, g1, g2;
mag_t C, D;
acb_init(m);
acb_init(t);
acb_init(log_t);
acb_init(u);
acb_init(v);
acb_init(g0);
acb_init(g1);
acb_init(g2);
mag_init(C);
mag_init(D);
/* m = mid(x) */
acb_set(m, x);
mag_zero(arb_radref(acb_realref(m)));
mag_zero(arb_radref(acb_imagref(m)));
/* t = alpha + im */
acb_mul_onei(t, m);
acb_add(t, t, alpha, prec);
/* log_t = log(alpha+im) */
acb_log(log_t, t, prec);
/* u = x - m */
acb_sub(u, x, m, prec);
acb_get_mag(D, u);
/* g0 = g(m) = (n+1) log(log(alpha+im)) - 2 pi m */
acb_log(g0, log_t, prec);
acb_mul_fmpz(g0, g0, n1, prec);
acb_const_pi(v, prec);
acb_mul(v, v, m, prec);
acb_mul_2exp_si(v, v, 1);
acb_sub(g0, g0, v, prec);
/* g1 = g'(m) (x-m); g'(m) = i(n+1)/((alpha+im) log(alpha+im)) - 2 pi */
acb_mul(g1, t, log_t, prec);
acb_inv(g1, g1, prec);
acb_mul_fmpz(g1, g1, n1, prec);
acb_mul_onei(g1, g1);
acb_const_pi(t, prec); /* recycle t */
acb_mul_2exp_si(t, t, 1);
acb_sub(g1, g1, t, prec);
acb_mul(g1, g1, u, prec);
/* absolute value -- optional (does not seem to affect speed much) */
acb_abs(acb_realref(g1), g1, prec);
arb_zero(acb_imagref(g1));
stieltjes_bound_quadratic_term(acb_realref(g2), x, n1, alpha, prec);
acb_exp(g0, g0, prec);
acb_exp(g1, g1, prec);
acb_exp(g2, g2, prec);
acb_get_mag(C, g0);
acb_get_mag(D, g1);
mag_mul(C, C, D);
acb_get_mag(D, g2);
mag_mul(C, C, D);
mag_mul_ui(C, C, 5);
arb_zero(B);
arf_set_mag(arb_midref(B), C);
acb_clear(m);
acb_clear(t);
acb_clear(log_t);
acb_clear(u);
acb_clear(v);
mag_clear(C);
mag_clear(D);
acb_clear(g0);
acb_clear(g1);
acb_clear(g2);
}
static void
stieltjes_bound_large(acb_t res, const acb_t x,
const fmpz_t n1, const acb_t alpha, slong prec)
{
arb_t B;
mag_t t;
arb_init(B);
mag_init(t);
prec = FLINT_MIN(prec, 30 + fmpz_bits(n1));
stieltjes_bound_large3(B, x, n1, alpha, prec);
arb_get_mag(t, B);
acb_zero(res);
arb_add_error_mag(acb_realref(res), t);
arb_add_error_mag(acb_imagref(res), t);
arb_clear(B);
mag_clear(t);
}
static void
stieltjes_integrand(acb_t res, const acb_t x, const fmpz_t n1,
const acb_t alpha, int analytic, slong prec)
{
acb_t t, u;
acb_init(t);
acb_init(u);
acb_mul_onei(t, x);
acb_add(t, t, alpha, prec);
acb_set_ui(u, 5);
/* check branch cut of logarithm */
if (arb_contains_zero(acb_imagref(t)) && !arb_is_positive(acb_realref(t)))
{
acb_indeterminate(res);
}
else if ((analytic || acb_rel_accuracy_bits(t) < prec - 10)
&& arb_gt(acb_realref(x), acb_realref(u)))
{
/* note: requires re(x) > 1 */
stieltjes_bound_large(res, x, n1, alpha, prec);
}
else
{
acb_const_pi(u, prec);
acb_mul(u, u, x, prec);
acb_sech(u, u, prec);
if (acb_is_finite(u))
{
acb_mul(u, u, u, prec);
acb_log(t, t, prec);
acb_pow_fmpz(t, t, n1, prec);
acb_mul(res, t, u, prec);
}
else
{
acb_indeterminate(res);
}
}
acb_clear(t);
acb_clear(u);
}
typedef struct
{
const fmpz * n1;
const acb_struct * alpha;
}
_stieltjes_param;
/* Compute the approximate saddle point. */
static void
stieltjes_omega(acb_t omega, const fmpz_t n1, const acb_t alpha, slong prec)
{
acb_t t, u;
fmpz_t k;
acb_init(t);
acb_init(u);
fmpz_init(k);
/* u = (n+1)i / (2pi) */
arb_set_fmpz(acb_imagref(t), n1);
acb_const_pi(u, prec);
acb_mul_2exp_si(u, u, 1);
acb_div(u, t, u, prec);
/* u / W(u) */
acb_lambertw(t, u, k, 0, prec);
acb_div(t, u, t, prec);
acb_sub(t, alpha, t, prec);
acb_mul_onei(t, t);
acb_set(omega, t);
acb_clear(t);
acb_clear(u);
fmpz_clear(k);
}
/* Compute an approximation of the magnitude of gamma_n. */
static void
stieltjes_mag_approx(arb_t C, mag_t tol, const fmpz_t n1, const acb_t alpha)
{
slong prec;
acb_t w, v, q;
prec = 32 + 2 * fmpz_bits(n1);
acb_init(w);
acb_init(v);
acb_init(q);
stieltjes_omega(w, n1, alpha, prec);
stieltjes_integrand(v, w, n1, alpha, 0, prec);
acb_set_fmpz(q, n1);
acb_sqrt(q, q, prec);
acb_mul(v, v, q, prec);
acb_get_mag(tol, v);
arb_set(C, acb_imagref(w));
mag_zero(arb_radref(C));
acb_clear(w);
acb_clear(v);
acb_clear(q);
}
int
_f_stieltjes(acb_ptr res, const acb_t x, void * param, slong order, slong prec)
{
const fmpz * n1;
const acb_struct * alpha;
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
n1 = ((const _stieltjes_param *) param)->n1;
alpha = ((const _stieltjes_param *) param)->alpha;
stieltjes_integrand(res, x, n1, alpha, order != 0, prec);
return 0;
}
/* Estimate log2(|gamma_n|), for fast cancellation estimate when n is small. */
static double
stieltjes_mag(double n)
{
double u, v, A, B;
double va, vb, t;
double pi = 3.141592653589793;
int i;
va = 1e-6;
vb = 0.5 * pi - 1e-6;
for (i = 0; i < 53; i++) /* bisection */
{
v = (va + vb) * 0.5;
t = 2 * pi * exp(v * tan(v)) - n * cos(v) / v;
if (t < 0.0)
va = v;
else
vb = v;
}
v = va;
u = v * tan(v);
A = 0.5 * log(u * u + v * v) - u / (u * u + v * v);
B = 2 * sqrt(2 * pi) * sqrt(u * u + v * v) * pow((u + 1) * (u + 1) + v * v, -0.25);
t = (log(B) + n * A - 0.5 * log(n)) / log(2);
return t;
}
static double _hypot(double x, double y) { return sqrt(x * x + y * y); }
/* log2 magnitude of integrand at z = x+yi; alpha = a+bi */
static double
integrand_mag(double n, double x, double y, double a, double b)
{
double t, u;
t = log(_hypot(a - y, b + x));
u = atan2(b + x, a - y);
t = log(_hypot(t,u)) * (n+1) - 2.0 * 3.1415926535897932 * x;
return t * 1.44269504088896341;
}
static double
find_x_maximizing_mag(double n, double y)
{
double xa, xb, xma, xmb, ma, mb;
int i;
xa = 1.0;
xb = n;
for (i = 0; i < 80; i++)
{
xma = xa + (xb - xa) / 3;
xmb = xb - (xb - xa) / 3;
ma = integrand_mag(n, xma, y, 0.5, 0.0);
mb = integrand_mag(n, xmb, y, 0.5, 0.0);
if (ma < mb)
xa = xma;
else
xb = xmb;
}
return xa;
}
static void
stieltjes_choose_N(arb_t N, const fmpz_t n1, const acb_t alpha, slong prec)
{
if (fmpz_bits(n1) < 30)
{
double nn, NN, aa, bb;
nn = fmpz_get_d(n1) - 1.0;
NN = FLINT_MAX(nn, 4);
aa = arf_get_d(arb_midref(acb_realref(alpha)), ARF_RND_DOWN);
bb = arf_get_d(arb_midref(acb_imagref(alpha)), ARF_RND_DOWN);
while (integrand_mag(nn, NN, 0.0, aa, bb) > -prec - 20)
{
NN *= 2.0;
if (NN > 1e30)
break;
}
arb_set_d(N, NN);
}
else
{
/* todo: account for huge alpha? */
arb_set_fmpz(N, n1);
}
}
static void
stieltjes_tail_bound(mag_t bound, const arb_t N, const fmpz_t n1, const acb_t alpha)
{
slong prec;
arb_t x, y, D;
acb_t aNi, logaNi;
mag_t t, u;
prec = MAG_BITS + fmpz_bits(n1);
arb_init(x);
arb_init(y);
arb_init(D);
acb_init(aNi);
acb_init(logaNi);
mag_init(t);
mag_init(u);
/* alpha + Ni */
acb_set(aNi, alpha);
arb_add(acb_imagref(aNi), acb_imagref(aNi), N, prec);
/* log(alpha + Ni) */
acb_log(logaNi, aNi, prec);
/* check (n+1)/(|alpha+Ni| log(alpha + Ni)|) < 2 */
acb_get_mag_lower(t, aNi);
acb_get_mag_lower(u, logaNi);
mag_mul_lower(t, t, u);
mag_inv(t, t);
mag_mul_fmpz(t, t, n1);
/* also check N >= |im(alpha)| + 2 */
arb_abs(x, acb_imagref(alpha));
arb_add_ui(x, x, 2, prec);
if (mag_cmp_2exp_si(t, 1) >= 0 || !arb_ge(N, x))
{
mag_inf(bound);
}
else
{
/* exp(-2 pi N) */
arb_set(x, N);
arb_mul_2exp_si(x, x, 1);
arb_const_pi(y, prec);
arb_mul(y, y, x, prec);
arb_neg(y, y);
arb_exp(y, y, prec);
/* |log|^(n+1) */
acb_get_mag(t, logaNi);
arf_set_mag(arb_midref(x), t);
mag_zero(arb_radref(x));
arb_pow_fmpz(x, x, n1, prec);
arb_mul(x, x, y, prec);
arb_get_mag(bound, x);
}
arb_clear(x);
arb_clear(y);
arb_clear(D);
acb_clear(aNi);
acb_clear(logaNi);
mag_clear(t);
mag_clear(u);
}
void
_acb_dirichlet_stieltjes_integral2(acb_t res, const fmpz_t n, const acb_t alpha, slong prec)
{
double gamma_mag, max_mag, cancellation, xa;
fmpz_t n1;
acb_t a, b, v, w;
slong wp;
mag_t tol, bound;
acb_calc_integrate_opt_t opt;
_stieltjes_param param;
/* integration points */
arb_t M, N, C;
/* required for the integral representation to be valid */
if (!arb_is_positive(acb_realref(alpha)))
{
acb_indeterminate(res);
return;
}
fmpz_init(n1);
arb_init(M);
arb_init(N);
arb_init(C);
acb_init(a);
acb_init(b);
acb_init(v);
acb_init(w);
mag_init(tol);
mag_init(bound);
fmpz_add_ui(n1, n, 1);
printf("stieltjes_integral2 "); fmpz_print(n); printf(" "); acb_printn(alpha, 30, 0); printf("%ld\n", prec);
param.n1 = n1;
param.alpha = alpha;
arb_set_ui(M, 10);
stieltjes_choose_N(N, n1, alpha, prec);
stieltjes_tail_bound(bound, N, n1, alpha);
printf("N = "); arb_printn(N, 30, 0); printf("\n");
printf("bound = "); mag_printd(bound, 10); printf("\n");
if (acb_is_real(alpha) &&
arf_cmpabs_2exp_si(arb_midref(acb_realref(alpha)), 2) < 0 &&
fmpz_cmp_ui(n1, 5000) < 0)
{
double nn = fmpz_get_ui(n1) - 1;
gamma_mag = stieltjes_mag(nn);
xa = find_x_maximizing_mag(nn, 0.0);
max_mag = integrand_mag(nn, xa, 0.0, 0.5, 0.0);
cancellation = FLINT_MAX(0, max_mag - gamma_mag);
if (cancellation < 10 + 0.1 * prec)
{
arb_zero(C);
mag_one(tol);
mag_mul_2exp_si(tol, tol, gamma_mag);
}
else
{
stieltjes_mag_approx(C, tol, n1, alpha);
cancellation = 0;
}
}
else
{
stieltjes_mag_approx(C, tol, n1, alpha);
cancellation = 0;
}
printf("C = "); arb_printn(C, 30, 0); printf("\n");
printf("cancellation = %f\n", cancellation);
mag_mul_2exp_si(tol, tol, -prec - 5);
/* todo: 1 * fmpz_bits(n1) should be enough, but acb powering is
inaccurate */
wp = prec + 2 * fmpz_bits(n1) + cancellation + 10;
printf("wp = %ld\n", wp);
acb_calc_integrate_opt_init(opt);
/* opt->verbose = 1; */
opt->deg_limit = 100 + 1.2 * prec; /* Small speedup. */
opt->verbose = 2;
if (arb_is_zero(C))
{
acb_zero(a);
acb_set_arb(b, N);
acb_calc_integrate(w, _f_stieltjes, &param, a, b, wp, tol, opt, wp);
acb_add(v, v, w, wp);
}
else
{
acb_zero(a); /* a = 0 */
acb_set_arb(b, M); /* b = M */
acb_calc_integrate(w, _f_stieltjes, &param, a, b, wp, tol, opt, wp);
acb_add(v, v, w, wp);
acb_set(a, b);
acb_set_arb(b, M);
arb_set(acb_imagref(b), C); /* b = M + i C */
acb_calc_integrate(w, _f_stieltjes, &param, a, b, wp, tol, opt, wp);
acb_add(v, v, w, wp);
acb_set(a, b);
arb_set(acb_realref(b), N); /* b = N + i C */
acb_calc_integrate(w, _f_stieltjes, &param, a, b, wp, tol, opt, wp);
acb_add(v, v, w, wp);
acb_set(a, b);
arb_zero(acb_imagref(b)); /* b = N */
acb_calc_integrate(w, _f_stieltjes, &param, a, b, wp, tol, opt, wp);
acb_add(v, v, w, wp);
}
printf("v = "); acb_printn(v, 30, 0); printf("\n");
acb_add_error_mag(v, bound);
acb_const_pi(b, wp);
acb_mul(v, v, b, wp);
acb_div_fmpz(v, v, n1, wp);
acb_neg(v, v);
if (acb_is_real(alpha))
arb_zero(acb_imagref(v));
acb_set_round(res, v, prec);
fmpz_clear(n1);
acb_clear(a);
acb_clear(b);
acb_clear(v);
acb_clear(w);
mag_clear(tol);
mag_clear(bound);
arb_clear(M);
arb_clear(N);
arb_clear(C);
}
void
_acb_dirichlet_stieltjes_integral(acb_t res, const fmpz_t n, const acb_t a, slong prec)
{
acb_t alpha;
acb_init(alpha);
/* alpha = a-1/2 */
acb_set_d(alpha, 0.5);
acb_sub(alpha, a, alpha, prec);
if (acb_is_real(a))
{
acb_conj(alpha, alpha);
_acb_dirichlet_stieltjes_integral2(res, n, alpha, prec);
}
else
{
acb_t r1, r2;
acb_init(r1);
acb_init(r2);
_acb_dirichlet_stieltjes_integral2(r1, n, alpha, prec);
acb_conj(alpha, alpha);
_acb_dirichlet_stieltjes_integral2(r2, n, alpha, prec);
acb_conj(r2, r2);
acb_add(res, r1, r2, prec);
acb_mul_2exp_si(res, res, -1);
acb_clear(r1);
acb_clear(r2);
}
acb_clear(alpha);
}
void
acb_dirichlet_stieltjes_integral(acb_t res, const fmpz_t n, const acb_t a, slong prec)
{
/* the algorithm only works for re(a) > 1/2 */
if (arf_cmp_si(arb_midref(acb_realref(a)), 1) < 0)
{
slong k, m, wp;
acb_t ak, t, s;
if (arf_cmp_si(arb_midref(acb_realref(a)), -prec) < 0)
{
acb_indeterminate(res);
return;
}
m = 1 - arf_get_si(arb_midref(acb_realref(a)), ARF_RND_FLOOR);
acb_init(ak);
acb_init(t);
acb_init(s);
wp = prec + 2 * fmpz_bits(n);
for (k = 0; k < m; k++)
{
acb_add_si(ak, a, k, wp);
acb_log(t, ak, wp);
acb_pow_fmpz(t, t, n, wp);
acb_div(t, t, ak, wp);
acb_add(s, s, t, wp);
}
acb_add_si(ak, a, m, wp);
_acb_dirichlet_stieltjes_integral(t, n, ak, prec);
acb_add(res, s, t, prec);
acb_clear(s);
acb_clear(t);
acb_clear(ak);
}
else
{
_acb_dirichlet_stieltjes_integral(res, n, a, prec);
}
}
void
acb_dirichlet_stieltjes_em(acb_t res, const fmpz_t n, const acb_t a, slong prec)
{
if (fmpz_cmp_ui(n, 10000) > 0)
{
acb_indeterminate(res);
}
else
{
slong nn, wp;
acb_ptr z;
acb_t s;
nn = fmpz_get_ui(n);
acb_init(s);
z = _acb_vec_init(nn + 1);
wp = prec * 1.05 + 2.2*nn + 10; /* todo: bogus */
/* todo: we don't want to compute all the coefficients */
acb_one(s);
_acb_poly_zeta_cpx_series(z, s, a, 1, nn + 1, wp);
arb_fac_ui(acb_realref(s), nn, prec + 10);
acb_mul_arb(res, z + nn, acb_realref(s), prec);
if (fmpz_is_odd(n))
acb_neg(res, res);
acb_clear(s);
_acb_vec_clear(z, nn + 1);
}
}
void
acb_dirichlet_stieltjes(acb_t res, const fmpz_t n, const acb_t a, slong prec)
{
slong cutoff;
if (acb_is_one(a) && fmpz_is_zero(n))
{
arb_const_euler(acb_realref(res), prec);
arb_zero(acb_imagref(res));
return;
}
if (fmpz_sgn(n) < 0)
{
flint_printf("stieltjes constants only defined for n >= 0");
flint_abort();
}
/* undefined at a = 0, -1, -2, ... */
if (acb_contains_int(a) && !arb_is_positive(acb_realref(a)))
{
acb_indeterminate(res);
return;
}
cutoff = FLINT_MAX(100, prec / 2);
cutoff = FLINT_MIN(cutoff, 10000);
if (fmpz_cmp_ui(n, cutoff) >= 0)
{
acb_dirichlet_stieltjes_integral(res, n, a, prec);
}
else
{
acb_dirichlet_stieltjes_em(res, n, a, prec);
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment