Created
September 13, 2018 05:09
-
-
Save fredrik-johansson/0a1edb681144b0de977984e059ef4fb9 to your computer and use it in GitHub Desktop.
debug stieltjes.c
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/* | |
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, ¶m, 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, ¶m, 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, ¶m, 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, ¶m, 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, ¶m, 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