Skip to content

Commit 1abaf57

Browse files
fredrik-johanssonfredrik-johansson
committed
first implementation of Lerch transcendent
1 parent b401d7c commit 1abaf57

File tree

6 files changed

+1120
-0
lines changed

6 files changed

+1120
-0
lines changed

‎acb_dirichlet.h‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,10 @@ void acb_dirichlet_zeta_jet(acb_t res, const acb_t s, int deflate, slong len, sl
4848

4949
void acb_dirichlet_hurwitz(acb_t res, const acb_t s, const acb_t a, slong prec);
5050

51+
void acb_dirichlet_lerch_phi_integral(acb_t res, const acb_t z, const acb_t s, const acb_t a, slong prec);
52+
void acb_dirichlet_lerch_phi_direct(acb_t res, const acb_t z, const acb_t s, const acb_t a, slong prec);
53+
void acb_dirichlet_lerch_phi(acb_t res, const acb_t z, const acb_t s, const acb_t a, slong prec);
54+
5155
void acb_dirichlet_stieltjes(acb_t res, const fmpz_t n, const acb_t a, slong prec);
5256

5357
typedef struct

‎acb_dirichlet/lerch_phi.c‎

Lines changed: 153 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,153 @@
1+
/*
2+
Copyright (C) 2022 Fredrik Johansson
3+
4+
This file is part of Arb.
5+
6+
Arb is free software: you can redistribute it and/or modify it under
7+
the terms of the GNU Lesser General Public License (LGPL) as published
8+
by the Free Software Foundation; either version 2.1 of the License, or
9+
(at your option) any later version. See <http://www.gnu.org/licenses/>.
10+
*/
11+
12+
#include "acb_hypgeom.h"
13+
#include "acb_dirichlet.h"
14+
15+
void
16+
acb_dirichlet_lerch_phi(acb_t res, const acb_t z, const acb_t s, const acb_t a, slong prec)
17+
{
18+
if (!acb_is_finite(z) || !acb_is_finite(s) || !acb_is_finite(a))
19+
{
20+
acb_indeterminate(res);
21+
return;
22+
}
23+
24+
if (acb_contains_int(a) && !arb_is_positive(acb_realref(a)))
25+
{
26+
if (!(acb_is_int(s) && arb_is_nonpositive(acb_realref(s))))
27+
{
28+
acb_indeterminate(res);
29+
return;
30+
}
31+
}
32+
33+
if (acb_is_zero(z))
34+
{
35+
acb_t t;
36+
acb_init(t);
37+
acb_neg(t, s);
38+
acb_pow(res, a, t, prec);
39+
acb_clear(t);
40+
return;
41+
}
42+
43+
if (acb_is_one(z))
44+
{
45+
arb_t one;
46+
arb_init(one);
47+
if (arb_gt(acb_realref(s), one))
48+
acb_dirichlet_hurwitz(res, s, a, prec);
49+
else
50+
acb_indeterminate(res);
51+
arb_clear(one);
52+
return;
53+
}
54+
55+
if (acb_equal_si(z, -1))
56+
{
57+
if (acb_is_one(a))
58+
{
59+
acb_dirichlet_eta(res, s, prec);
60+
}
61+
else if (acb_is_one(s))
62+
{
63+
/* (psi((a+1)/2) - psi(a/2))/2 */
64+
acb_t t, u;
65+
acb_init(t);
66+
acb_init(u);
67+
acb_mul_2exp_si(t, a, -1);
68+
acb_digamma(t, t, prec);
69+
acb_add_ui(u, a, 1, prec);
70+
acb_mul_2exp_si(u, u, -1);
71+
acb_digamma(u, u, prec);
72+
acb_sub(res, u, t, prec);
73+
acb_mul_2exp_si(res, res, -1);
74+
acb_clear(t);
75+
acb_clear(u);
76+
}
77+
else
78+
{
79+
/* 2^(-s) (zeta(s,a/2) - zeta(s,(a+1)/2)) */
80+
acb_t t, u;
81+
acb_init(t);
82+
acb_init(u);
83+
acb_mul_2exp_si(t, a, -1);
84+
acb_hurwitz_zeta(t, s, t, prec);
85+
acb_add_ui(u, a, 1, prec);
86+
acb_mul_2exp_si(u, u, -1);
87+
acb_hurwitz_zeta(u, s, u, prec);
88+
acb_sub(t, t, u, prec);
89+
acb_neg(u, s);
90+
acb_set_ui(res, 2);
91+
acb_pow(res, res, u, prec);
92+
acb_mul(res, res, t, prec);
93+
acb_clear(t);
94+
acb_clear(u);
95+
}
96+
return;
97+
}
98+
99+
if (acb_is_zero(s))
100+
{
101+
acb_sub_ui(res, z, 1, prec + 5);
102+
acb_neg(res, res);
103+
acb_inv(res, res, prec);
104+
return;
105+
}
106+
107+
if (acb_is_one(s))
108+
{
109+
acb_t t, u;
110+
acb_init(t);
111+
acb_init(u);
112+
acb_one(t);
113+
acb_add_ui(u, a, 1, prec + 5);
114+
acb_hypgeom_2f1(t, t, a, u, z, ACB_HYPGEOM_2F1_BC, prec + 5);
115+
acb_div(res, t, a, prec);
116+
if (!acb_is_finite(res))
117+
acb_indeterminate(res);
118+
acb_clear(t);
119+
acb_clear(u);
120+
return;
121+
}
122+
123+
if (acb_is_one(a) && !acb_contains_zero(z))
124+
{
125+
acb_t t;
126+
acb_init(t);
127+
acb_polylog(t, s, z, prec);
128+
acb_div(res, t, z, prec);
129+
acb_clear(t);
130+
return;
131+
}
132+
133+
{
134+
mag_t zm, lim;
135+
mag_init(zm);
136+
mag_init(lim);
137+
138+
acb_get_mag(zm, z);
139+
mag_set_d(lim, 0.75);
140+
141+
if (mag_cmp(zm, lim) <= 0)
142+
{
143+
acb_dirichlet_lerch_phi_direct(res, z, s, a, prec);
144+
}
145+
else
146+
{
147+
acb_dirichlet_lerch_phi_integral(res, z, s, a, prec);
148+
}
149+
150+
mag_clear(zm);
151+
mag_clear(lim);
152+
}
153+
}

‎acb_dirichlet/lerch_phi_direct.c‎

Lines changed: 150 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,150 @@
1+
/*
2+
Copyright (C) 2022 Fredrik Johansson
3+
4+
This file is part of Arb.
5+
6+
Arb is free software: you can redistribute it and/or modify it under
7+
the terms of the GNU Lesser General Public License (LGPL) as published
8+
by the Free Software Foundation; either version 2.1 of the License, or
9+
(at your option) any later version. See <http://www.gnu.org/licenses/>.
10+
*/
11+
12+
#include "acb_dirichlet.h"
13+
14+
void
15+
acb_dirichlet_lerch_phi_direct(acb_t res, const acb_t z, const acb_t s, const acb_t a, slong prec)
16+
{
17+
slong N, Nmax, wp, n;
18+
int a_real;
19+
acb_t negs, t, u, sum;
20+
mag_t C, S, zmag, tail_bound, tm, tol;
21+
22+
if (!acb_is_finite(z) || !acb_is_finite(s) || !acb_is_finite(a))
23+
{
24+
acb_indeterminate(res);
25+
return;
26+
}
27+
28+
if (acb_contains_int(a) && !arb_is_positive(acb_realref(a)))
29+
{
30+
if (!(acb_is_int(s) && arb_is_nonpositive(acb_realref(s))))
31+
{
32+
acb_indeterminate(res);
33+
return;
34+
}
35+
}
36+
37+
acb_init(negs);
38+
acb_init(t);
39+
acb_init(u);
40+
acb_init(sum);
41+
42+
acb_neg(negs, s);
43+
44+
mag_init(C);
45+
mag_init(S);
46+
mag_init(zmag);
47+
mag_init(tail_bound);
48+
mag_init(tm);
49+
mag_init(tol);
50+
51+
a_real = acb_is_real(a);
52+
wp = prec + 10;
53+
54+
acb_get_mag(zmag, z);
55+
56+
/* first term: 1/a^s */
57+
acb_pow(sum, a, negs, wp);
58+
59+
acb_get_mag(tol, sum);
60+
mag_mul_2exp_si(tol, tol, -wp);
61+
62+
if (a_real)
63+
{
64+
/* Tail bound |z|^N / |(a+N)^s| * sum C^k, C = |z| * exp(max(0, -re(s)) / (a+N)) */
65+
arb_nonnegative_part(acb_realref(t), acb_realref(negs));
66+
arb_get_mag(S, acb_realref(t));
67+
}
68+
else
69+
{
70+
/* Tail bound |z|^N / |(a+N)^s| * sum C^k, C = |z| * exp(|s / (a+N)|) */
71+
acb_get_mag(S, s);
72+
}
73+
74+
Nmax = 100 * prec + 0.1 * prec * n_sqrt(prec);
75+
Nmax = FLINT_MAX(Nmax, 1);
76+
Nmax = FLINT_MIN(Nmax, WORD_MAX / 2);
77+
78+
mag_inf(tail_bound);
79+
80+
for (N = 1; N <= Nmax; N = FLINT_MAX(N+4, N*1.1))
81+
{
82+
acb_add_ui(t, a, N, 53);
83+
84+
if (arb_is_positive(acb_realref(t)))
85+
{
86+
acb_get_mag_lower(C, t);
87+
mag_div(C, S, C);
88+
mag_exp(C, C);
89+
mag_mul(C, C, zmag);
90+
mag_geom_series(C, C, 0);
91+
92+
if (mag_is_finite(C))
93+
{
94+
mag_pow_ui(tail_bound, zmag, N);
95+
mag_mul(tail_bound, tail_bound, C);
96+
acb_pow(t, t, negs, 53);
97+
acb_get_mag(C, t);
98+
mag_mul(tail_bound, tail_bound, C);
99+
100+
if (mag_cmp(tail_bound, tol) <= 0)
101+
break;
102+
}
103+
else
104+
{
105+
mag_inf(tail_bound);
106+
}
107+
}
108+
}
109+
110+
if (mag_is_finite(tail_bound))
111+
{
112+
acb_one(t);
113+
114+
for (n = 1; n < N; n++)
115+
{
116+
if (n % 8 == 0 && !acb_is_real(z))
117+
acb_pow_ui(t, z, n, wp);
118+
else
119+
acb_mul(t, t, z, wp);
120+
121+
acb_add_ui(u, a, n, wp);
122+
acb_pow(u, u, negs, wp);
123+
acb_mul(u, t, u, wp);
124+
acb_add(sum, sum, u, wp);
125+
}
126+
127+
if (acb_is_real(z) && acb_is_real(s) && acb_is_real(a))
128+
arb_add_error_mag(acb_realref(sum), tail_bound);
129+
else
130+
acb_add_error_mag(sum, tail_bound);
131+
132+
acb_set_round(res, sum, prec);
133+
}
134+
else
135+
{
136+
acb_indeterminate(res);
137+
}
138+
139+
mag_clear(C);
140+
mag_clear(S);
141+
mag_clear(zmag);
142+
mag_clear(tail_bound);
143+
mag_clear(tm);
144+
mag_clear(tol);
145+
146+
acb_clear(negs);
147+
acb_clear(t);
148+
acb_clear(u);
149+
acb_clear(sum);
150+
}

0 commit comments

Comments
 (0)