This section describes routines related to L-functions. We first introduce
the basic concept and notations, then explain how to represent them in GP.
Let Γ_{ℝ}(s) = π^{-s/2}Γ(s/2), where Γ is Euler's gamma
function. Given d ≥ 1 and a d-tuple A = [α_{1},...,α_{d}] of
complex numbers, we let γ_{A}(s) = ∏_{α ∈ A} Γ_{ℝ}(s +
α).

Given a sequence a = (a_{n})_{n ≥ 1} of complex numbers (such that a_{1} = 1),
a positive *conductor* N ∈ ℤ, and a *gamma factor*
γ_{A} as above, we consider the Dirichlet series
L(a,s) = ∑_{n ≥ 1} a_{n} n^{-s}
and the attached completed function
Λ(a,s) = N^{s/2}γ_{A}(s).L(a,s).

Such a datum defines an *L-function* if it satisfies the three
following assumptions:

***** [Convergence] The a_{n} = O_ε(n^{k1+ε}) have polynomial
growth, equivalently L(s) converges absolutely in some right half-plane
Re(s) > k_{1} + 1.

***** [Analytic continuation] L(s) has a meromorphic continuation to the
whole complex plane with finitely many poles.

***** [Functional equation] There exist an integer k, a complex number
ε (usually of modulus 1), and an attached sequence a^{*}
defining both an L-function L(a^{*},s) satisfying the above two assumptions
and a completed function Λ(a^{*},s) = N^{s/2}γ_{A}(s).
L(a^{*},s), such that
Λ(a,k-s) = ε Λ(a^{*},s)
for all regular points.

More often than not in number theory we have a^ *= a (which
forces |ε |= 1), but this needs not be the case. If a is a real
sequence and a = a^{*}, we say that L is *self-dual*. We do not assume
that the a_{n} are multiplicative, nor equivalently that L(s) has an Euler
product.

**Remark.**
Of course, a determines the L-function, but the (redundant) datum a,a^{*},
A, N, k, ε describes the situation in a form more suitable for fast
computations; knowing the polar part r of Λ(s) (a rational function
such that Λ-r is holomorphic) is also useful. A subset of these,
including only finitely many a_{n}-values will still completely determine L
(in suitable families), and we provide routines to try and compute missing
invariants from whatever information is available.

**Important Caveat.**
The implementation assumes that the implied constants in the O_ε are
small. In our generic framework, it is impossible to return proven results
without more detailed information about the L function. The intended use of
the L-function package is not to prove theorems, but to experiment and
formulate conjectures, so all numerical results should be taken with a grain
of salt. One can always increase `realbitprecision`

and recompute: the
difference estimates the actual absolute error in the original output.

**Note.** The requested precision has a major impact on runtimes.
Because of this, most L-function routines, in particular `lfun`

itself,
specify the requested precision in *bits*, not in decimal digits.
This is transparent for the user once `realprecision`

or
`realbitprecision`

are set. We advise to manipulate precision via
`realbitprecision`

as it allows finer granularity: `realprecision`

increases by increments of 64 bits, i.e. 19 decimal digits at a time.

Given an L-function as above, we define an attached theta function
via Mellin inversion: for any positive real t > 0, we let
θ(a,t) := (1)/(2π i)∫_{Re(s) = c} t^{-s} Λ(s) ds
where c is any positive real number c > k_{1}+1 such that c + Re(a) > 0
for all a ∈ A. In fact, we have
θ(a,t) = ∑_{n ≥ 1} a_{n} K(nt/N^{1/2})
where
K(t) := (1)/(2π i)∫_{Re(s) = c} t^{-s} γ_{A}(s) ds.
Note that this function is analytic and actually makes sense for complex t,
such that Re(t^{2/d}) > 0, i.e. in a cone containing the positive real
half-line. The functional equation for Λ translates into
θ(a,1/t) - ε t^kθ(a^{*},t) = P_Λ(t),
where P_Λ is an explicit polynomial in t and log t given by the
Taylor development of the polar part of Λ: there are no log's if
all poles are simple, and P = 0 if Λ is entire. The values
θ(t) are generally easier to compute than the L(s), and this
functional equation provides a fast way to guess possible values for
missing invariants in the L-function definition.

We have 3 levels of description:

***** an `Lmath`

is an arbitrary description of the underlying
mathematical situation (to which e.g., we associate the a_{p} as traces of
Frobenius elements); this is done via constructors to be described in the
subsections below.

***** an `Ldata`

is a computational description of situation, containing
the complete datum (a,a^{*},A,k,N,ε,r). Where a and a^{*} describe
the coefficients (given n,b we must be able to compute [a_{1},...,a_{n}]
with bit accuracy b), A describes the Euler factor, the (classical) weight
is k, N is the conductor, and r describes the polar part of L(s).
This is obtained via the function `lfuncreate`

. N.B. For motivic
L-functions, the motivic weight w is w = k-1; but we also support
non-motivic L-functions.

**Design problem.** All components of an `Ldata`

should be given
exactly since the accuracy to which they must be computed is not bounded a
priori; but this is not always possible, in particular for ε and r.

***** an `Linit`

contains an `Ldata`

and everything needed for fast
*numerical* computations. It specifies the functions to be considered
(either L^{(j)}(s) or θ^{(j)}(t) for derivatives of order j ≤
m, where m is now fixed) and specifies a *domain* which limits
the range of arguments (t or s, respectively to certain cones and
rectangular regions) and the output accuracy. This is obtained via the
functions `lfuninit`

or `lfunthetainit`

.

All the functions which are specific to L or theta functions share the
prefix `lfun`

. They take as first argument either an `Lmath`

, an
`Ldata`

, or an `Linit`

. If a single value is to be computed,
this makes no difference, but when many values are needed (e.g. for plots or
when searching for zeros), one should first construct an `Linit`

attached to the search range and use it in all subsequent calls.
If you attempt to use an `Linit`

outside the range for which it was
initialized, a warning is issued, because the initialization is
performed again, a major inefficiency:

? Z = lfuncreate(1); \\ Riemann zeta ? L = lfuninit(Z, [1/2, 0, 100]); \\ zeta(1/2+it), |t| < 100 ? lfun(L, 1/2) \\ OK, within domain %3 = -1.4603545088095868128894991525152980125 ? lfun(L, 0) \\ not on critical line ! *** lfun: Warning: lfuninit: insufficient initialization. %4 = -0.50000000000000000000000000000000000000 ? lfun(L, 1/2, 1) \\ attempt first derivative ! *** lfun: Warning: lfuninit: insufficient initialization. %5 = -3.9226461392091517274715314467145995137

For many L-functions, passing from `Lmath`

to an `Ldata`

is
inexpensive: in that case one may use `lfuninit`

directly from the
`Lmath`

even when evaluations in different domains are needed. The
above example could equally have skipped the `lfuncreate`

:

? L = lfuninit(1, [1/2, 0, 100]); \\ zeta(1/2+it), |t| < 100

In fact, when computing a single value, you can even skip
`lfuninit`

:

? L = lfun(1, 1/2, 1); \\ zeta'(1/2) ? L = lfun(1, 1+x+O(x^5)); \\ first 5 terms of Taylor development at 1

Both give the desired results with no warning.

**Complexity.**
The implementation requires O(N(|t|+1))^{1/2} coefficients a_{n}
to evaluate L of conductor N at s = σ + i t.

We now describe the available high-level constructors, for built-in L functions.

Given a Dirichlet character χ:(ℤ/Nℤ)^{*} → ℂ, we let
L(χ, s) = ∑_{n ≥ 1} χ(n) n^{-s}.
Only primitive characters are supported. Given a fundamental discriminant
D, the function L((D/.), s), for the quadratic Kronecker symbol, is encoded
by the `t_INT`

D. This includes Riemann ζ function via the special
case D = 1.

More general characters can be represented in a variety of ways:

***** via Conrey notation (see `znconreychar`

): χ_{N}(m,.)
is given as the `t_INTMOD`

`Mod(m,N)`

.

***** via a *znstar* structure describing the abelian group (ℤ/Nℤ)^{*},
where the character is given in terms of the *znstar* generators:

? G = znstar(100, 1); \\ (Z/100Z)^{*}? G.cyc \\ ~ Z/20 . g1 + Z/2 . g2 for some generators g1 and g2 %2 = [20, 2] ? G.gen %3 = [77, 51] ? chi = [a, b] \\ maps g1 to e(a/20) and g2 to e(b/2); e(x) = exp(2ipi x)

More generally, let (ℤ/Nℤ)^ *= ⨁ (ℤ/d_{i}ℤ) g_{i} be given via a
*znstar* structure G (`G.cyc`

gives the d_{i} and `G.gen`

the
g_{i}). A *character* χ on G is given by a row vector
v = [a_{1},...,a_{n}] such that χ(∏ g_{i}^{ni}) = exp(2π i∑ a_{i}
n_{i} / d_{i}). The pair [G, v] encodes the *primitive* character
attached to χ.

***** in fact, this construction [G, m] describing a character
is more general: m is also allowed to be a Conrey index as seen above,
or a Conrey logarithm (see `znconreylog`

), and the latter format is
actually the fastest one.

***** it is also possible to view Dirichlet characters as Hecke characters
over K = ℚ (see below), for a modulus [N, [1]] but this is both more
complicated and less efficient.

In all cases, a non-primitive character is replaced by the attached primitive character.

The Dedekind zeta function of a number field K = ℚ[X]/(T) is encoded
either by the defining polynomial T, or any absolute number fields
structure (preferably at least a *bnf*).

Given a finite order Hecke character χ: Cl_{f}(K) → ℂ, we let
L(χ, s) = ∑_{A ⊂ OK} χ(A) (N_{K/ℚ}A)^{-s}.

Let Cl_{f}(K) = ⨁ (ℤ/d_{i}ℤ) g_{i} given by a *bnr* structure with
generators: the d_{i} are given by `K.cyc`

and the g_{i} by `K.gen`

.
A *character* χ on the ray class group is given by a row vector
v = [a_{1},...,a_{n}] such that χ(∏ g_{i}^{ni}) = exp(2π i∑
a_{i} n_{i} / d_{i}). The pair [*bnr*, v] encodes the *primitive*
character attached to χ.

? K = bnfinit(x^2-60); ? Cf = bnrinit(K, [7, [1,1]], 1); \\ f = 7 oo_{1}oo_{2}? Cf.cyc %3 = [6, 2, 2] ? Cf.gen %4 = [[2, 1; 0, 1], [22, 9; 0, 1], [-6, 7]~] ? lfuncreate([Cf, [1,0,0]]); \\ χ(g_{1}) = ζ_{6}, χ(g_{2}) = χ(g_{3}) = 1

Dirichlet characters on (ℤ/Nℤ)^{*} are a special case,
where K = ℚ:

? Q = bnfinit(x); ? Cf = bnrinit(Q, [100, [1]]); \\ for odd characters on (Z/100Z)*

For even characters, replace by `bnrinit(K, N)`

. Note that the simpler
direct construction in the previous section will be more efficient.

Given a Galois number field N/ℚ with group G = `galoisinit`

(N),
a representation ρ of G over the cyclotomic field ℚ(ζ_{n})
is specified by the matrices giving the images of `G.gen`

by ρ.
The corresponding Artin L function is created using `lfunartin`

.

P = quadhilbert(-47); \\ degree 5, Galois group D_{5}N = nfinit(nfsplitting(P)); \\ Galois closure G = galoisinit(N); [s,t] = G.gen; \\ order 5 and 2 L = lfunartin(N,G, [[a,0;0,a^-1],[0,1;1,0]], 5); \\ irr. degree 2

In the above, the polynomial variable (here `a`

) represents
ζ_{5} := exp(2iπ/5) and the two matrices give the images of
s and t. Here, priority of `a`

must be lower than the priority
of `x`

.

L-function of elliptic curves over number fields are supported.

? E = ellinit([1,1]); ? L = lfuncreate(E); \\ L-function of E/Q ? E2 = ellinit([1,a], nfinit(a^2-2)); ? L2 = lfuncreate(E2); \\ L-function of E/Q(sqrt(2))

L-function of hyperelliptic genus-2 curve can be created with
`lfungenus2`

. To create the L function of the curve
y^2+(x^3+x^2+1)y = x^2+x:

? L = lfungenus2([x^2+x, x^3+x^2+1]);

Currently, the model needs to be minimal at 2, and if the conductor is even, its valuation at 2 might be incorrect (a warning is issued).

An eta quotient is created by applying `lfunetaquo`

to a matrix with
2 columns [m, r_{m}] representing
f(τ) := ∏_{m} η(mτ)^{rm}.
It is currently assumed that f is a self-dual cuspidal form on
Γ_{0}(N) for some N.
For instance, the L-function ∑ τ(n) n^{-s}
attached to Ramanujan's Δ function is encoded as follows

? L = lfunetaquo(Mat([1,24])); ? lfunan(L, 100) \\ first 100 values of tau(n)

More general modular forms defined by modular symbols will be added later.

When no direct constructor is available, you can still input an L function
directly by supplying [a, a^{*},A, k, N, ε, r] to `lfuncreate`

(see `??lfuncreate`

for details).

It is *strongly* suggested to first check consistency of the created
L-function:

? L = lfuncreate([a, as, A, k, N, eps, r]); ? lfuncheckfeq(L) \\ check functional equation

Compute the L-function value L(s), or if `D`

is set, the
derivative of order `D`

at s. The parameter
`L`

is either an Lmath, an Ldata (created by `lfuncreate`

, or an
Linit (created by `lfuninit`

), preferrably the latter if many values
are to be computed.

The argument s is also allowed to be a power series; for instance, if s =
α + x + O(x^n), the function returns the Taylor expansion of order n
around α. The result is given with absolute error less than 2^{-B},
where B = realbitprecision.

**Caveat.** The requested precision has a major impact on runtimes.
It is advised to manipulate precision via `realbitprecision`

as
explained above instead of `realprecision`

as the latter allows less
granularity: `realprecision`

increases by increments of 64 bits, i.e. 19
decimal digits at a time.

? lfun(x^2+1, 2) \\ Lmath: Dedekind zeta for Q(i) at 2 %1 = 1.5067030099229850308865650481820713960 ? L = lfuncreate(ellinit("5077a1")); \\ Ldata: Hasse-Weil zeta function ? lfun(L, 1+x+O(x^4)) \\ zero of order 3 at the central point %3 = 0.E-58 - 5.[...] E-40*x + 9.[...] E-40*x^2 + 1.7318[...]*x^3 + O(x^4) \\ Linit: zeta(1/2+it), |t| < 100, and derivative ? L = lfuninit(1, [100], 1); ? T = lfunzeros(L, [1,25]); %5 = [14.134725[...], 21.022039[...]] ? z = 1/2 + I*T[1]; ? abs( lfun(L, z) ) %7 = 8.7066865533412207420780392991125136196 E-39 ? abs( lfun(L, z, 1) ) %8 = 0.79316043335650611601389756527435211412 \\ simple zero

The library syntax is `GEN `

.**lfun0**(GEN L, GEN s, long D, long bitprec)

Returns the `Linit`

structure attached to the Dedekind zeta function
of the number field L (see `lfuninit`

), given a subfield K such that
L/K is abelian.
Here `polrel`

defines L over K, as usual with the priority of the
variable of `bnfK`

lower than that of `polrel`

.
`sdom`

and `der`

are as in `lfuninit`

.

? D = -47; K = bnfinit(y^2-D); ? rel = quadhilbert(D); T = rnfequation(K.pol, rel); \\ degree 10 ? L = lfunabelianrelinit(T,K,rel, [2,0,0]); \\ at 2 time = 84 ms. ? lfun(L, 2) %4 = 1.0154213394402443929880666894468182650 ? lfun(T, 2) \\ using parisize > 300MB time = 652 ms. %5 = 1.0154213394402443929880666894468182656

As the example shows, using the (abelian) relative structure is more efficient than a direct computation. The difference becomes drastic as the absolute degree increases while the subfield degree remains constant.

The library syntax is `GEN `

.**lfunabelianrelinit**(GEN bnfL, GEN bnfK, GEN polrel, GEN sdom, long der, long bitprec)

Compute the first n terms of the Dirichlet series attached to the
L-function given by `L`

(`Lmath`

, `Ldata`

or `Linit`

).

? lfunan(1, 10) \\ Riemann zeta %1 = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1] ? lfunan(5, 10) \\ Dirichlet L-function for kronecker(5,.) %2 = [1, -1, -1, 1, 0, 1, -1, -1, 1, 0]

The library syntax is `GEN `

.**lfunan**(GEN L, long n, long prec)

Returns the `Ldata`

structure attached to the
Artin L-function provided by the representation ρ of the Galois group
of the extension K/ℚ, defined over the cyclotomic field ℚ(ζ_{n}),
where *nf* is the nfinit structure attached to K,
*gal* is the galoisinit structure attached to K/ℚ, and *rho* is
given either

***** by the values of its character on the conjugacy classes
(see `galoisconjclasses`

and `galoischartable`

)

***** or by the matrices that are the images of the generators

.*gal*.gen

Cyclotomic numbers in `rho`

are represented by polynomials, whose
variable is understood as the complex number exp(2 i π/n).

In the following example we build the Artin L-functions attached to the two
irreducible degree 2 representations of the dihedral group D_{10} defined
over ℚ(ζ_{5}), for the extension H/ℚ where H is the Hilbert class
field of ℚ(sqrt{-47}).
We show numerically some identities involving Dedekind ζ functions and
Hecke L series.

? P = quadhilbert(-47) %1 = x^5 + 2*x^4 + 2*x^3 + x^2 - 1 ? N = nfinit(nfsplitting(P)); ? G = galoisinit(N); \\ D_10 ? [T,n] = galoischartable(G); ? T \\ columns give the irreducible characters %5 = [1 1 2 2] [1 -1 0 0] [1 1 -y^3 - y^2 - 1 y^3 + y^2] [1 1 y^3 + y^2 -y^3 - y^2 - 1] ? n %6 = 5 ? L2 = lfunartin(N,G, T[,2], n); ? L3 = lfunartin(N,G, T[,3], n); ? L4 = lfunartin(N,G, T[,4], n); ? s = 1 + x + O(x^4); ? lfun(-47,s) - lfun(L2,s) %11 ~ 0 ? lfun(1,s)*lfun(-47,s)*lfun(L3,s)^2*lfun(L4,s)^2 - lfun(N,s) %12 ~ 0 ? lfun(1,s)*lfun(L3,s)*lfun(L4,s) - lfun(P,s) %13 ~ 0 ? bnr = bnrinit(bnfinit(x^2+47),1,1); ? bnr.cyc %15 = [5] \\ Z/5Z: 4 non-trivial ray class characters ? lfun([bnr,[1]], s) - lfun(L3, s) %16 ~ 0 ? lfun([bnr,[2]], s) - lfun(L4, s) %17 ~ 0 ? lfun([bnr,[3]], s) - lfun(L3, s) %18 ~ 0 ? lfun([bnr,[4]], s) - lfun(L4, s) %19 ~ 0

The first identity identifies the non-trivial abelian character with
(-47,.); the second is the factorization of the regular representation of
D_{10}; the third is the factorization of the natural representation of
D_{10} ⊂ S_{5}; and the final four are the expressions of the degree 2
representations as induced from degree 1 representations.

The library syntax is `GEN `

.**lfunartin**(GEN nf, GEN gal, GEN rho, long n, long bitprec)

Given the data attached to an L-function (`Lmath`

, `Ldata`

or `Linit`

), check whether the functional equation is satisfied.
This is most useful for an `Ldata`

constructed "by hand", via
`lfuncreate`

, to detect mistakes.

If the function has poles, the polar part must be specified. The routine
returns a bit accuracy b such that |w - ^{w}| < 2^{b}, where w is
the root number contained in `data`

, and ^{w} is a computed value
derived from θ(t) and θ(1/t) at some t != 0 and
the assumed functional equation. Of course, a large negative value of the
order of `realbitprecision`

is expected.

If t is given, it should be close to the unit disc for efficiency and such that θ(t) != 0. We then check the functional equation at that t.

? \pb 128 \\ 128 bits of accuracy ? default(realbitprecision) %1 = 128 ? L = lfuncreate(1); \\ Riemann zeta ? lfuncheckfeq(L) %3 = -124

i.e. the given data is consistent to within 4 bits for the particular check consisting of estimating the root number from all other given quantities. Checking away from the unit disc will either fail with a precision error, or give disappointing results (if θ(1/t) is large it will be computed with a large absolute error)

? lfuncheckfeq(L, 2+I) %4 = -115 ? lfuncheckfeq(L,10) *** at top-level: lfuncheckfeq(L,10) *** ^ — — — — — — *** lfuncheckfeq: precision too low in lfuncheckfeq.

The library syntax is `long `

.**lfuncheckfeq**(GEN L, GEN t = NULL, long bitprec)

Compute the conductor of the given L-function
(if the structure contains a conductor, it is ignored);
`ab`

= [a,b] is the interval where we expect to find the conductor;
it may be given as a single scalar b, in which case we look in [1,b].
Increasing `ab`

slows down the program but gives better accuracy for the
result.

If `flag`

is 0 (default), give either the conductor found as an
integer, or a vector (possibly empty) of conductors found. If `flag`

is
1, same but give the computed floating point approximations to the
conductors found, without rounding to integers. It `flag`

is 2, give
all the conductors found, even those far from integers.

**Caveat.** This is a heuristic program and the result is not
proven in any way:

? L = lfuncreate(857); \\ Dirichlet L function for kronecker(857,.) ? \p19 realprecision = 19 significant digits ? lfunconductor(L) %2 = [17, 857] ? lfunconductor(L,,1) \\ don't round %3 = [16.99999999999999999, 857.0000000000000000] ? \p38 realprecision = 38 significant digits ? lfunconductor(L) %4 = 857

**Note.** This program should only be used if the primes dividing the
conductor are unknown, which is rare. If they are known, a direct
search through possible prime exponents using `lfuncheckfeq`

will
be more efficient and rigorous:

? E = ellinit([0,0,0,4,0]); /* Elliptic curve y^2 = x^3+4x */ ? E.disc \\ |disc E| = 2^12 %2 = -4096 \\ create Ldata by hand. Guess that root number is 1 and conductor N ? L(N) = lfuncreate([n->ellan(E,n), 0, [0,1], 2, N, 1]); ? fordiv(E.disc, d, print(d,": ",lfuncheckfeq(L(d)))) 1: 0 2: 0 4: -1 8: -2 16: -3 32: -127 64: -3 128: -2 256: -2 512: -1 1024: -1 2048: 0 4096: 0 ? lfunconductor(L(1)) \\ lfunconductor ignores conductor = 1 in Ldata ! %5 = 32

The above code assumed that root number was 1;
had we set it to -1, none of the `lfuncheckfeq`

values would have been
acceptable:

? L2(N) = lfuncreate([n->ellan(E,n), 0, [0,1], 2, N, -1]); ? [ lfuncheckfeq(L2(d)) | d<-divisors(E.disc) ] %7 = [0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, -1, -1]

The library syntax is `GEN `

.**lfunconductor**(GEN L, GEN ab = NULL, long 10000], long bitprec)

Estimate the cost of running
`lfuninit(L,sdom,der)`

at current bit precision. Returns [t,b], to
indicate that t coefficients a_{n} will be computed, as well as t values of
`gammamellininv`

, all at bit accuracy b.
A subsequent call to `lfun`

at s evaluates a polynomial of degree t
at exp(h s) for some real parameter h, at the same bit accuracy b.
If L is already an `Linit`

, then *sdom* and *der* are ignored
and are best left omitted; the bit accuracy is also inferred from L: in
short we get an estimate of the cost of using that particular `Linit`

.

? \pb 128 ? lfuncost(1, [100]) \\ for zeta(1/2+I*t), |t| < 100 %1 = [7, 242] \\ 7 coefficients, 242 bits ? lfuncost(1, [1/2, 100]) \\ for zeta(s) in the critical strip, |Im s| < 100 %2 = [7, 246] \\ now 246 bits ? lfuncost(1, [100], 10) \\ for zeta(1/2+I*t), |t| < 100 %3 = [8, 263] \\ 10th derivative increases the cost by a small amount ? lfuncost(1, [10^5]) %3 = [158, 113438] \\ larger imaginary part: huge accuracy increase ? L = lfuncreate(polcyclo(5)); \\ Dedekind zeta for Q(zeta_{5}) ? lfuncost(L, [100]) \\ at s = 1/2+I*t), |t| < 100 %5 = [11457, 582] ? lfuncost(L, [200]) \\ twice higher %6 = [36294, 1035] ? lfuncost(L, [10^4]) \\ much higher: very costly ! %7 = [70256473, 45452] ? \pb 256 ? lfuncost(L, [100]); \\ doubling bit accuracy %8 = [17080, 710]

In fact, some L functions can be factorized algebraically
by the `lfuninit`

call, e.g. the Dedekind zeta function of abelian
fields, leading to much faster evaluations than the above upper bounds.
In that case, the function returns a vector of costs as above for each
individual function in the product actually evaluated:

? L = lfuncreate(polcyclo(5)); \\ Dedekind zeta for Q(zeta_{5}) ? lfuncost(L, [100]) \\ a priori cost %2 = [11457, 582] ? L = lfuninit(L, [100]); \\ actually perform all initializations ? lfuncost(L) %4 = [[16, 242], [16, 242], [7, 242]]

The Dedekind function of this abelian quartic field is the product of four Dirichlet L-functions attached to the trivial character, a non-trivial real character and two complex conjugate characters. The non-trivial characters happen to have the same conductor (hence same evaluation costs), and correspond to two evaluations only since the two conjugate characters are evaluated simultaneously. For a total of three L-functions evaluations, which explains the three components above. Note that the actual cost is much lower than the a priori cost in this case.

The library syntax is `GEN `

.
Also available is
**lfuncost0**(GEN L, GEN sdom = NULL, long der, long bitprec)`GEN `

when L is **lfuncost**(GEN L, GEN dom, long der, long bitprec)*not* an `Linit`

; the return value is a `t_VECSMALL`

in this case.

This low-level routine creates `Ldata`

structures, needed by
*lfun* functions, describing an L-function and its functional equation.
You are urged to use a high-level constructor when one is available,
and this function accepts them, see `??lfun`

:

? L = lfuncreate(1); \\ Riemann zeta ? L = lfuncreate(5); \\ Dirichlet L-function for quadratic character (5/.) ? L = lfuncreate(x^2+1); \\ Dedekind zeta for Q(i) ? L = lfuncreate(ellinit([0,1])); \\ L-function of E/Q: y^2=x^3+1

One can then use, e.g., `Lfun(L,s)`

to directly
evaluate the respective L-functions at s, or `lfuninit(L, [c,w,h]`

to initialize computations in the rectangular box Re(s-c) ≤ w,
Im(s) ≤ h.

We now describe the low-level interface, used to input non-builtin L-functions. The input is now a 6 or 7 component vector V = [a, astar, Vga, k, N, eps, poles], whose components are as follows:

***** `V[1] = a`

encodes the Dirichlet series coefficients (a_{n}). The
preferred format is a closure of arity 1: `n- > vector(n,i,a(i))`

giving
the vector of the first n coefficients. The closure is allowed to return
a vector of more than n coefficients (only the first n will be
considered) or even less than n, in which case loss of accuracy will occur
and a warning that `#an`

is less than expected is issued. This
allows to precompute and store a fixed large number of Dirichlet
coefficients in a vector v and use the closure `n- > v`

, which
does not depend on n. As a shorthand for this latter case, you can input
the vector v itself instead of the closure.

? z = lfuncreate([n->vector(n,i,1), 1, [0], 1, 1, 1, 1]); \\ Riemann zeta ? lfun(z,2) - Pi^2/6 %2 = -5.877471754111437540 E-39

A second format is limited to L-functions affording an
Euler product. It is a closure of arity 2 `(p,d)- > F(p)`

giving the
local factor L_{p}(X) at p as a rational function, to be evaluated at
p^{-s} as in `direuler`

; d is set to `logint`

(n,p) + 1, where
n is the total number of Dirichlet coefficients (a_{1},...,a_{n}) that will
be computed. In other words, the smallest integer d such that p^d > n.
This parameter d allows to compute only part of
L_{p} when p is large and L_{p} expensive to compute: any polynomial
(or `t_SER`

) congruent to L_{p} modulo X^d is acceptable since only
the coefficients of X^0,..., X^{d-1} are needed to expand the Dirichlet
series. The closure can of course ignore this parameter:

? z = lfuncreate([(p,d)->1/(1-x), 1, [0], 1, 1, 1, 1]); \\ Riemann zeta ? lfun(z,2) - Pi^2/6 %4 = -5.877471754111437540 E-39

One can describe separately the generic local factors coefficients
and the bad local factors by setting `dir`

= [F, L_{bad}],
were L_{bad} = [[p_{1},L_{p1}],...,[p_{k},L_{pk}]], where F
describes the generic local factors as above, except that when p = p_{i}
for some i ≤ k, the coefficient a_{p} is directly set to L_{pi}
instead of calling F.

N = 15; E = ellinit([1, 1, 1, -10, -10]); \\ = "15a1" F(p,d) = 1 / (1 - ellap(E,p)*'x + p*'x^2); Lbad = [[3, 1/(1+'x)], [5, 1/(1-'x)]]; L = lfuncreate([[F,Lbad], 0, [0,1], 2, N, ellrootno(E)]);

Of course, in this case, `lfuncreate(E)`

is preferable!

***** `V[2] = astar`

is the Dirichlet series coefficients of the dual
function, encoded as `a`

above. The sentinel values 0 and 1 may
be used for the special cases where a = a^{*} and a = a^{*},
respectively.

***** `V[3] = Vga`

is the vector of α_{j} such that the gamma
factor of the L-function is equal to
γ_{A}(s) = ∏_{1 ≤ j ≤ d}Γ_{ℝ}(s+α_{j}),
where Γ_{ℝ}(s) = π^{-s/2}Γ(s/2).
This same syntax is used in the `gammamellininv`

functions.
In particular the length d of `Vga`

is the degree of the L-function.
In the present implementation, the α_{j} are assumed to be exact
rational numbers. However when calling theta functions with *complex*
(as opposed to real) arguments, determination problems occur which may
give wrong results when the α_{j} are not integral.

***** `V[4] = k`

is a positive integer k. The functional equation relates
values at s and k-s. For instance, for an Artin L-series such as a
Dedekind zeta function we have k = 1, for an elliptic curve k = 2, and
for a modular form, k is its weight. For motivic L-functions, the
*motivic* weight w is w = k-1.

By default we assume that a_{n} = O_ε(n^{k1+ε}), where
k_{1} = w and even k_{1} = w/2 when the L function has no pole
(Ramanujan-Petersson). If this is not the case, you can replace the
k argument by a vector [k,k_{1}], where k_{1} is the upper bound you can
assume.

***** `V[5] = N`

is the conductor, an integer N ≥ 1, such that
Λ(s) = N^{s/2}γ_{A}(s)L(s) with γ_{A}(s) as above.

***** `V[6] = eps`

is the root number ϵ, i.e., the
complex number (usually of modulus 1) such that
Λ(a, k-s) = ϵ Λ(a^{*}, s).

***** The last optional component `V[7] = poles`

encodes the poles of the
L or Λ-functions, and is omitted if they have no poles.
A polar part is given by a list of 2-component vectors
[β,P_{β}(x)], where
β is a pole and the power series P_{β}(x) describes
the attached polar part, such that L(s) - P_β(s-β) is holomorphic
in a neighbourhood of β. For instance P_β = r/x+O(1) for a
simple pole at β or r_{1}/x^2+r_{2}/x+O(1) for a double pole.
The type of the list describing the polar part allows to distinguish between
L and Λ: a `t_VEC`

is attached to L, and a `t_COL`

is attached to Λ. Unless a = a^{*} (coded by `astar`

equal to 0 or 1), it is mandatory to specify the polar part of Λ
rather than those of L since the poles of L^{*} cannot be infered from the
latter ! Whereas the functional equation allows to deduce the polar part of
Λ^{*} from the polar part of Λ.

Finally, if a = a^{*}, we allow a shortcut to describe
the frequent situation where L has at most simple pole, at s = k,
with residue r a complex scalar: you may then input `poles`

= r.
This value r can be set to 0 if unknown and it will be computed.

The library syntax is `GEN `

.**lfuncreate**(GEN obj)

Creates the `Ldata`

structure (without initialization) corresponding
to the quotient of the Dirichlet series L_{1} and L_{2} given by
`L1`

and `L2`

. Assume that v_{z}(L_{1}) ≥ v_{z}(L_{2}) at all
complex numbers z: the construction may not create new poles, nor increase
the order of existing ones.

The library syntax is `GEN `

.**lfundiv**(GEN L1, GEN L2, long bitprec)

Returns the `Ldata`

structure attached to the L function
attached to the modular form
z` ⟼ `

∏_{i = 1}^n η(M_{i,1} z)^{Mi,2}
It is currently assumed that f is a self-dual cuspidal form on
Γ_{0}(N) for some N.
For instance, the L-function ∑ τ(n) n^{-s}
attached to Ramanujan's Δ function is encoded as follows

? L = lfunetaquo(Mat([1,24])); ? lfunan(L, 100) \\ first 100 values of tau(n)

For convenience, a `t_VEC`

is also accepted instead of
a factorization matrix with a single row:

? L = lfunetaquo([1,24]); \\ same as above

The library syntax is `GEN `

.**lfunetaquo**(GEN M)

Returns the `Ldata`

structure attached to the L function
attached to the genus-2 curve defined by y^2 = F(x) or
y^2+Q(x) y = P(x) if F = [P,Q].
Currently, the model needs to be minimal at 2, and if the conductor
is even, its valuation at 2 might be incorrect (a warning is issued).

The library syntax is `GEN `

.**lfungenus2**(GEN F)

Variant of the Hardy Z-function given by `L`

, used for
plotting or locating zeros of L(k/2+it) on the critical line.
The precise definition is as
follows: if as usual k/2 is the center of the critical strip, d is the
degree, α_{j} the entries of `Vga`

giving the gamma factors,
and ϵ the root number, then if we set
s = k/2+it = ρ e^{iθ} and
E = (d(k/2-1)+∑_{1 ≤ j ≤ d}α_{j})/2, the computed function at t is
equal to
Z(t) = ϵ^{-1/2}Λ(s).|s|^{-E}e^{dtθ/2} ,
which is a real function of t for self-dual Λ,
vanishing exactly when L(k/2+it) does on the critical line. The
normalizing factor |s|^{-E}e^{dtθ/2} compensates the
exponential decrease of γ_{A}(s) as t → oo so that
Z(t) ~ 1.

? T = 100; \\ maximal height ? L = lfuninit(1, [T]); \\ initialize for zeta(1/2+it), |t|<T ? \p19 \\ no need for large accuracy ? ploth(t = 0, T, lfunhardy(L,t))

Using `lfuninit`

is critical for this particular
applications since thousands of values are computed. Make sure to initialize
up to the maximal t needed: otherwise expect to see many warnings for
unsufficient initialization and suffer major slowdowns.

The library syntax is `GEN `

.**lfunhardy**(GEN L, GEN t, long bitprec)

Initalization function for all functions linked to the
computation of the L-function L(s) encoded by `L`

, where
s belongs to the rectangular domain `sdom`

= [*center*,w,h]
centered on the real axis, |Re(s)-*center*| ≤ w, |Im(s)| ≤ h,
where all three components of `sdom`

are real and w, h are
non-negative. `der`

is the maximum order of derivation that will be used.
The subdomain [k/2, 0, h] on the critical line (up to height h)
can be encoded as [h] for brevity. The subdomain [k/2, w, h]
centered on the critical line can be encoded as [w, h] for brevity.

The argument `L`

is an `Lmath`

, an `Ldata`

or an `Linit`

. See
`??Ldata`

and `??lfuncreate`

for how to create it.

The height h of the domain is a *crucial* parameter: if you only
need L(s) for real s, set h to 0.
The running time is roughly proportional to
(B / d+π h/4)^{d/2+3}N^{1/2},
where B is the default bit accuracy, d is the degree of the
L-function, and N is the conductor (the exponent d/2+3 is reduced
to d/2+2 when d = 1 and d = 2). There is also a dependency on w,
which is less crucial, but make sure to use the smallest rectangular
domain that you need.

? L0 = lfuncreate(1); \\ Riemann zeta ? L = lfuninit(L0, [1/2, 0, 100]); \\ for zeta(1/2+it), |t| < 100 ? lfun(L, 1/2 + I) ? L = lfuninit(L0, [100]); \\ same as above !

The library syntax is `GEN `

.**lfuninit0**(GEN L, GEN sdom, long der, long bitprec)

Compute the completed L-function Λ(s) = N^{s/2}γ(s)L(s),
or if `D`

is set, the derivative of order `D`

at s.
The parameter `L`

is either an `Lmath`

, an `Ldata`

(created by
`lfuncreate`

, or an `Linit`

(created by `lfuninit`

), preferrably the
latter if many values are to be computed.

The result is given with absolute error less than 2^{-B}|γ(s)N^{s/2}|,
where B = realbitprecision.

The library syntax is `GEN `

.**lfunlambda0**(GEN L, GEN s, long D, long bitprec)

Returns `[valeven,valodd,omminus,omplus]`

,
where `valeven`

(resp., `valodd`

) is the vector of even (resp., odd)
periods of the modular form given by `L`

, and `omminus`

and
`omplus`

the corresponding real numbers ω^- and ω^+
normalized in a noncanonical way. For the moment, only for modular forms of even weight.

The library syntax is `GEN `

.**lfunmfspec**(GEN L, long bitprec)

Creates the `Ldata`

structure (without initialization) corresponding
to the product of the Dirichlet series given by `L1`

and
`L2`

.

The library syntax is `GEN `

.**lfunmul**(GEN L1, GEN L2, long bitprec)

Computes the order of the possible zero of the L-function at the center k/2 of the critical strip; return 0 if L(k/2) does not vanish.

If m is given and has a non-negative value, assumes the order is at most m. Otherwise, the algorithm chooses a sensible default:

***** if the L argument is an `Linit`

, assume that a multiple zero at
s = k / 2 has order less than or equal to the maximal allowed derivation
order.

***** else assume the order is less than 4.

You may explicitly increase this value using optional argument m; this
overrides the default value above. (Possibly forcing a recomputation
of the `Linit`

.)

The library syntax is `long `

.**lfunorderzero**(GEN L, long m, long bitprec)

Returns the `Ldata`

structure attached to the Θ function
of the lattice attached to the definite positive quadratic form Q.

? L = lfunqf(matid(2)); ? lfunqf(L,2) %2 = 6.0268120396919401235462601927282855839 ? lfun(x^2+1,2)*4 %3 = 6.0268120396919401235462601927282855839

The library syntax is `GEN `

.**lfunqf**(GEN Q, long prec)

Given the `Ldata`

attached to an L-function (or the output of
`lfunthetainit`

), compute the root number and the residues.

The output is a 3-component vector
[[[a_{1},r_{1}],...,[a_{n}, r_{n}], [[b_{1}, R_{1}],...,[b_{m}, R_{m}]], w],
where r_{i} is the
polar part of L(s) at a_{i}, R_{i} is is the polar part of Λ(s) at
b_{i} or [0,0,r] if there is no pole,
and w is the root number. In the present implementation,

***** either the polar part must be completely known (and is then arbitrary):
the function determines the root number,

? L = lfunmul(1,1); \\ zeta^2 ? [r,R,w] = lfunrootres(L); ? r \\ single pole at 1, double %3 = [[1, 1.[...]*x^-2 + 1.1544[...]*x^-1 + O(x^0)]] ? w %4 = 1 ? R \\ double pole at 0 and 1 %5 = [[1,[...]], [0,[...]]

***** or at most a single pole is allowed: the function computes both
the root number and the residue (0 if no pole).

The library syntax is `GEN `

.**lfunrootres**(GEN data, long bitprec)

Returns the `Ldata`

structure attached to the L function
attached to the m-th symmetric power of the elliptic curve E defined over
the rationals.

The library syntax is `GEN `

.**lfunsympow**(GEN E, ulong m)

Compute the value of the m-th derivative
at t of the theta function attached to the L-function given by `data`

.
`data`

can be either the standard L-function data, or the output of
`lfunthetainit`

. The result is given with absolute error less than
2^{-B}, where B = realbitprecision.

The theta function is defined by the formula
Θ(t) = ∑_{n ≥ 1}a(n)K(nt/sqrt(N)), where a(n) are the coefficients
of the Dirichlet series, N is the conductor, and K is the inverse Mellin
transform of the gamma product defined by the `Vga`

component.
Its Mellin transform is equal to Λ(s)-P(s), where Λ(s)
is the completed L-function and the rational function P(s) its polar part.
In particular, if the L-function is the L-function of a modular form
f(τ) = ∑_{n ≥ 0}a(n)q^n with q = exp(2π iτ), we have
Θ(t) = 2(f(it/sqrt{N})-a(0)). Note that a(0) = -L(f,0) in this case.

The library syntax is `GEN `

.**lfuntheta**(GEN data, GEN t, long m, long bitprec)

This function estimates the cost of running
`lfunthetainit(L,tdom,m)`

at current bit precision. Returns the number of
coefficients a_{n} that would be computed. This also estimates the
cost of a subsequent evaluation `lfuntheta`

, which must compute
that many values of `gammamellininv`

at the current bit precision.
If L is already an `Linit`

, then *tdom* and m are ignored
and are best left omitted: we get an estimate of the cost of using that
particular `Linit`

.

? \pb 1000 ? L = lfuncreate(1); \\ Riemann zeta ? lfunthetacost(L); \\ cost for theta(t), t real >= 1 %1 = 15 ? lfunthetacost(L, 1 + I); \\ cost for theta(1+I). Domain error ! *** at top-level: lfunthetacost(1,1+I) *** ^ — — — — — — -- *** lfunthetacost: domain error in lfunthetaneed: arg t > 0.785 ? lfunthetacost(L, 1 + I/2) \\ for theta(1+I/2). %2 = 23 ? lfunthetacost(L, 1 + I/2, 10) \\ for theta^((10))(1+I/2). %3 = 24 ? lfunthetacost(L, [2, 1/10]) \\ cost for theta(t), |t| >= 2, |arg(t)| < 1/10 %4 = 8 ? L = lfuncreate( ellinit([1,1]) ); ? lfunthetacost(L) \\ for t >= 1 %6 = 2471

The library syntax is `long `

.**lfunthetacost0**(GEN L, GEN tdom = NULL, long m, long bitprec)

Initalization function for evaluating the m-th derivative of theta
functions with argument t in domain *tdom*. By default (*tdom*
omitted), t is real, t ≥ 1. Otherwise, *tdom* may be

***** a positive real scalar ρ: t is real, t ≥ ρ.

***** a non-real complex number: compute at this particular t; this
allows to compute θ(z) for any complex z satisfying |z| ≥ |t|
and |arg z| ≤ |arg t|; we must have |2 arg z / d| < π/2, where
d is the degree of the Γ factor.

***** a pair [ρ,α]: assume that |t| ≥ ρ and |arg t| ≤
α; we must have |2α / d| < π/2, where d is the degree of
the Γ factor.

? \p500 ? L = lfuncreate(1); \\ Riemann zeta ? t = 1+I/2; ? lfuntheta(L, t); \\ direct computation time = 30 ms. ? T = lfunthetainit(L, 1+I/2); time = 30 ms. ? lfuntheta(T, t); \\ instantaneous

The T structure would allow to quickly compute θ(z) for any z in the cone delimited by t as explained above. On the other hand

? lfuntheta(T,I) *** at top-level: lfuntheta(T,I) *** ^ — — — — -- *** lfuntheta: domain error in lfunthetaneed: arg t > 0.785398163397448

The initialization is equivalent to

? lfunthetainit(L, [abs(t), arg(t)])

The library syntax is `GEN `

.**lfunthetainit**(GEN L, GEN tdom = NULL, long m, long bitprec)

Creates the Ldata structure (without initialization) corresponding to the
twist of L by the primitive character attached to the Dirichlet character
`chi`

. The conductor of the character must be coprime to the conductor
of the L-function L.

The library syntax is `GEN `

.**lfuntwist**(GEN L, GEN chi)

`lim`

being either a positive upper limit or a non-empty real
interval inside [0,+ oo [, computes an
ordered list of zeros of L(s) on the critical line up to the given
upper limit or in the given interval. Use a naive algorithm which may miss
some zeros: it assumes that two consecutive zeros at height T ≥ 1
differ at least by 2π/ω, where
ω := `divz`

.(dlog(T/2π) +d+ 2log(N/(π/2)^d)).
To use a finer search mesh, set divz to some integral value
larger than the default ( = 8).

? lfunzeros(1, 30) \\ zeros of Rieman zeta up to height 30 %1 = [14.134[...], 21.022[...], 25.010[...]] ? #lfunzeros(1, [100,110]) \\ count zeros with 100 <= Im(s) <= 110 %2 = 4

The algorithm also assumes that all zeros are simple except possibly on the real axis at s = k/2 and that there are no poles in the search interval. (The possible zero at s = k/2 is repeated according to its multiplicity.)

Should you pass an `Linit`

argument to the function, beware that the
algorithm needs at least

L = lfuninit(Ldata, [T+1])

where T is the upper bound of the interval defined by
`lim`

: this allows to detect zeros near T. Make sure that your
`Linit`

domain contains this one, i.e. a domain [1,T+1] is fine but
[0, T] is not! The algorithm assumes that a multiple zero at s = k / 2
has order less than or equal to the maximal derivation order allowed by the
`Linit`

. You may increase that value in the `Linit`

but this is
costly: only do it for zeros of low height or in `lfunorderzero`

instead.

The library syntax is `GEN `

.**lfunzeros**(GEN L, GEN lim, long divz, long bitprec)