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
nonmotivic L-functions.

**Technical note.** When some components of an `Ldata`

cannot be
given exactly, usually r or ε, the `Ldata`

may be given as a
*closure*. When evaluated at a given precision, the closure must return
all components as exact data or floating point numbers at the requested
precision, see `??lfuncreate`

. The reason for this technicality is that
the accuracy to which we must compute is not bounded a priori and unknown
at this stage: it depends on the domain where we evaluate the L-function.

***** 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 nonzero
integer D, the `t_INT`

D encodes the function L((D_{0}/.), s), for the
quadratic Kronecker symbol attached to the fundamental discriminant
D_{0} = `coredisc`

(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_{j}ℤ) g_{j} be given via a
*znstar* structure G (`G.cyc`

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

the
g_{j}). A *character* χ on G is given by a row vector
v = [a_{1},...,a_{n}] such that χ(∏_{j} g_{j}^{nj})
= exp(2π i∑_{j} a_{j} n_{j} / d_{j}). 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 label as seen above,
or a Conrey logarithm (see `znconreylog`

), and the latter format is
actually the fastest one. Instead
of a single character as above, one may use the construction
`lfuncreate([G, vchi])`

where `vchi`

is a nonempty vector of
characters *of the same conductor* (more precisely, whose attached
primitive characters have the same conductor) and *same parity*.
The function is then vector-valued, where the values are ordered as the
characters in `vchi`

. Conrey labels cannot be used in this last format
because of the need to distinguish a single character given by a row vector
of integers and a vector of characters given by their labels: use
`znconreylog(G,i)`

first to convert a label to Conrey logarithm.

***** 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 nonprimitive 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 (a *nf* is enough).

An alternative description for the Dedekind zeta function of an Abelian
extension of a number field is to use class-field theory parameters
[*bnr*, *subg*], see `bnrinit`

.

? bnf = bnfinit(a^2 - a - 9); ? bnr = bnrinit(bnf, [2, [0,0]]); subg = Mat(3); ? L = lfuncreate([bnr, subg]);

Let K be a number field given as a `bnfinit`

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

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

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

.
A *character* χ on the ray class group is given by a row vector
v = [a_{1},...,a_{n}] such that χ(∏_{j} g_{j}^{nj})
= exp(2π i∑_{j} a_{j} n_{j} / d_{j}). 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. Instead
of a single character as above, one may use the construction
`lfuncreate([Cf, vchi])`

where `vchi`

is a nonempty vector of
characters *of the same conductor* (more precisely, whose attached
primitive characters have the same conductor). The function is then
vector-valued, where the values are ordered as the characters in `vchi`

.

Given a Hecke *Grossencharacter* χ: \A^{×} → ℂ^{×} of
conductor 𝔣, we let
L(χ, s) = ∑_{A ⊂ ℤ_{K}, A+𝔣 = ℤ_{K}} χ(A) (N_{K/ℚ}A)^{-s}.

Let C_{K}(𝔪) = \A^{×}/(K^{×}.U(𝔪)) be an id\`ele class
group of modulus 𝔪 given by a *gchar* structure *gc* (see
`gcharinit`

and Section se:GCHAR).
A Grossencharacter χ on C_{K}(𝔪) is given by a row
vector of size `#`

.*gc*.cyc

? gc = gcharinit(x^3+4*x-1,[5,[1]]); \\ mod = 5.oo ? gc.cyc %3 = [4, 2, 0, 0] ? gcharlog(gc,idealprimedec(gc.bnf,5)[1]) \\ logarithm map C_{K}(𝔪) → ℝ^{n}? chi = [1,0,0,1,0]~; ? gcharduallog(gc,chi) \\ row vector of coefficients in ℝ^{n}? L = lfuncreate([gc,chi]); \\ non algebraic L-function ? lfunzeros(L,1) ? lfuneuler(L,2) \\ Euler factor at 2

Finite order Hecke characters are a special case.

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)

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 nontrivial 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 nontrivial 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} = θ(1/t) t^{-k} / θ(t) is a computed value
derived from the assumed functional equation.
Of course, the expected result is a large negative value of the order of
`realbitprecision`

. But if θ is very small
at t, you should first increase `realbitprecision`

by
-log_{2} |θ(t)|, which is
positive if θ is small, to get a meaningful result.
Note that t should be close to the unit disc for efficiency and such that
θ(t) ! = 0. If the parameter t is omitted, we check the
functional equation at the "random" complex number t = 335/339 + I/7.

? \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 case of Dedekind zeta functions.** Dedekind zeta function for
a number field K = ℚ[X]/(T) is in general computed
(assuming Artin conjecture) as (ζ_{K}/ζ_{k}) x ζ_{k},
where k is a
maximal subfield, applied recursively if possible. When K/ℚ is Galois,
the zeta function is directly decomposed as a product of Artin
L-functions.

These decompositions are computed when `lfuninit`

is called. The
behavior of `lfuncheckfeq`

is then different depending of its argument

***** the artificial query `lfuncheckfeq`

(T) serves little purpose
since we already know that the technical parameters are theoretically
correct; we just obtain an estimate on the accuracy they allow. This is
computed directly, without using the above decomposition. And is likely to
be very costly if the degree of T is large, possibly overflowing the
possibilities of the implementation.

***** a query `L = lfuninit(T, ...); lfuncheckfeq(L)`

on the other hand
returns the maximum of the `lfuncheckfeq`

values for all involved
L-functions, giving a general consistency check and again an estimate
for the accuracy of computed values.

At the default accuracy of 128 bits:

? T = polcyclo(43); ? lfuncheckfeq(T); *** at top-level: lfuncheckfeq(T) *** ^ — — — — — *** lfuncheckfeq: overflow in lfunthetacost. ? lfuncheckfeq(lfuninit(T, [2])) time = 107 ms. %2 = -122

The library syntax is `long `

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

Computes the conductor of the given L-function (if the structure
contains a conductor, it is ignored). Two methods are available,
depending on what we know about the conductor, encoded in the `setN`

parameter:

***** `setN`

is a scalar: we know nothing but expect that the conductor
lies in the interval [1, `setN`

].

If *flag* is 0 (default), gives either the conductor found as an
integer, or a vector (possibly empty) of conductors found. If *flag* is
1, same but gives the computed floating point approximations to the
conductors found, without rounding to integers. It *flag* is 2, gives
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

Increasing `setN`

or increasing `realbitprecision`

slows down the program but gives better accuracy for the result. This
algorithm should only be used if the primes dividing the conductor are
unknown, which is uncommon.

***** `setN`

is a vector of possible conductors; for instance
of the form `D1 * divisors(D2)`

, where D_{1} is the known part
of the conductor and D_{2} is a multiple of the contribution of the
bad primes.

In that case, *flag* is ignored and the routine uses `lfuncheckfeq`

.
It returns [N,e] where N is the best conductor in the list and e is the
value of `lfuncheckfeq`

for that N. When no suitable conductor exist or
there is a tie among best potential conductors, return the empty vector
`[]`

.

? 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]); \\ lfunconductor ignores conductor = 1 in Ldata ! ? lfunconductor(L(1), divisors(E.disc)) %5 = [32, -127] ? fordiv(E.disc, d, print(d,": ",lfuncheckfeq(L(d)))) \\ direct check 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

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 = lfuncreate([n->ellan(E,n), 0, [0,1], 2, 0, -1]); ? lfunconductor(L2, divisors(E.disc)) %7 = []

The library syntax is `GEN `

.**lfunconductor**(GEN L, GEN setN = NULL, long flag, long bitprec)

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

at current bit precision, given by a vector
[t, b].

***** If L contains the root number, 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 the root number is *not* known, then more values of a_{n} may
be needed in order to compute it, and the returned value of t takes this
into account (it may not be the exact value in this case but is always
an upper bound). Fewer than t `gammamellininv`

will be needed, and
a call to `lfun`

evaluates a polynomial of degree less that t, still
at 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`

.
Note that in this case, the root number is always already known and you get
the right value of t (corresponding to the number of past calls to
`gammamellinv`

and the actual degree of the evaluated polynomial).

? \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 is cheaper %8 = [17080, 710] ? \p38 ? K = bnfinit(y^2 - 4493); [P] = idealprimedec(K,1123); f = [P,[1,1]]; ? R = bnrinit(K, f); R.cyc %10 = [1122] ? L = lfuncreate([R, [7]]); \\ Hecke L-function ? L[6] %12 = 0 \\ unknown root number ? \pb 3000 ? lfuncost(L, [0], 1) %13 = [1171561, 3339] ? L = lfuninit(L, [0], 1); time = 1min, 56,426 ms. ? lfuncost(L) %14 = [826966, 3339]

In the final example, the root number was unknown and
extra coefficients a_{n} were needed to compute it (1171561). Once the
initialization is performed we obtain the lower value t = 826966, which
corresponds to the number of `gammamellinv`

computed and the actual
degree of the polynomial to be evaluated to compute a value within the
prescribed domain.

Finally, 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 nontrivial real character and two complex conjugate characters. The nontrivial 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.
We advise using a high-level constructor when one is available, see
`??lfun`

, and this function accepts them:

? 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 nonbuiltin 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_{p_{1}}],...,[p_{k},L_{p_{k}}]], 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_{p_{i}}
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.

**When one component is not exact.**
Alternatively, `obj`

can be a closure of arity 0 returning the above
vector to the current real precision. This is needed if some components
are not available exactly but only through floating point approximations.
The closure allows algorithms to recompute them to higher accuracy when
needed. Compare

? Ld1() = [n->lfunan(Mod(2,7),n),1,[0],1,7,((-13-3*sqrt(-3))/14)^(1/6)]; ? Ld2 = [n->lfunan(Mod(2,7),n),1,[0],1,7,((-13-3*sqrt(-3))/14)^(1/6)]; ? L1 = lfuncreate(Ld1); ? L2 = lfuncreate(Ld2); ? lfun(L1,1/2+I*200) \\ OK %5 = 0.55943925130316677665287870224047183265 - 0.42492662223174071305478563967365980756*I ? lfun(L2,1/2+I*200) \\ all accuracy lost %6 = 0.E-38 + 0.E-38*I

The accuracy lost in `Ld2`

is due to the root number being given to
an insufficient precision. To see what happens try

? Ld3() = printf("prec needed: %ld bits",getlocalbitprec());Ld1() ? L3 = lfuncreate(Ld3); prec needed: 64 bits ? z3 = lfun(L3,1/2+I*200) prec needed: 384 bits %16 = 0.55943925130316677665287870224047183265 - 0.42492662223174071305478563967365980756*I

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)

Creates the `Ldata`

structure (without initialization) corresponding
to the dual L-function ^{L} of L. If k and ϵ are
respectively the weight and root number of L, then the following formula
holds outside poles, up to numerical errors:
Λ(L, s) = ϵ Λ(^{L}, k - s).

? L = lfunqf(matdiagonal([1,2,3,4])); ? eps = lfunrootres(L)[3]; k = L[4]; ? M = lfundual(L); lfuncheckfeq(M) %3 = -127 ? s= 1+Pi*I; ? a = lfunlambda(L,s); ? b = eps * lfunlambda(M,k-s); ? exponent(a - b) %7 = -130

The library syntax is `GEN `

.**lfundual**(GEN L, 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)

Return the Euler factor at p of the
L-function given by `L`

(`Lmath`

, `Ldata`

or `Linit`

),
assuming the L-function admits an Euler product factorization and that it
can be determined.

? E=ellinit([1,3]); ? lfuneuler(E,7) %2 = 1/(7*x^2-2*x+1) ? L=lfunsympow(E,2); ? lfuneuler(L,11) %4 = 1/(-1331*x^3+275*x^2-25*x+1)

The library syntax is `GEN `

.**lfuneuler**(GEN L, GEN p, long prec)

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, 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: let k/2 be the center of the critical strip, d be the
degree, `Vga`

= (α_{j})_{j ≤ d} given the gamma factors,
and ϵ be the root number; we set
s = k/2+it = ρ e^{iθ} and
2E = d(k/2-1) + Re(∑_{1 ≤ j ≤ d}α_{j}). Assume first that
Λ is self-dual, then the computed function at t is equal to
Z(t) = ϵ^{-1/2}Λ(s).ρ^{-E}e^{dtθ/2} ,
which is a real function of t
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. For non-self-dual Λ, the definition is the same
except we drop the ϵ^{-1/2} term (which is not well defined since
it depends on the chosen dual sequence a^{*}(n)): Z(t) is still of the
order of 1 and still vanishes where L(k/2+it) does, but it needs no
longer be real-valued.

? 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
nonnegative. `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 !

**Riemann-Siegel formula.**
If L is a function of degree d = 1, then a completely different
algorithm is implemented which can compute with complexity N sqrt{h} (for
fixed accuracy B). So it handles larger imaginary parts than the default
implementation. But this variant is less efficient when the imaginary part of
s is tiny and the dependency in B is still in O(B^{2+1/2}).

For such functions, you can use *sdom* = `[]`

to indicate that you
are only interested in relatively high imaginary parts and do not want to
perform any initialization:

? L = lfuninit(1, []); \\ Riemann zeta ? #lfunzeros(L, [10^12, 10^12+1]) time = 1min, 31,496 ms. %2 = 4

If you ask instead for
`lfuninit(1, [10^12+1])`

, the initialization is restricted by some
cutoff value (depending on the conductor, but less than 10^4 in any case):
up to that point, the standard algorithm is used (and uses the
initialization); and above the cutoff, we switch to Riemann-Siegel. Note that
this is quite wasteful if only values with imaginary parts larger than 10^4
are needed.

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)

Let L be the L-function attached to a modular eigenform f of
weight k, as given by `lfunmf`

. In even weight, returns
`[ve,vo,om,op]`

, where `ve`

(resp., `vo`

) is the vector of even
(resp., odd) periods of f and `om`

and `op`

the corresponding
real numbers ω^{−} and ω^{+} normalized in a noncanonical way.
In odd weight `ominus`

is the same as `op`

and we
return `[v,op]`

where v is the vector of all periods.

? D = mfDelta(); mf = mfinit(D); L = lfunmf(mf, D); ? [ve, vo, om, op] = lfunmfspec(L) %2 = [[1, 25/48, 5/12, 25/48, 1], [1620/691, 1, 9/14, 9/14, 1, 1620/691],\ 0.0074154209298961305890064277459002287248,\ 0.0050835121083932868604942901374387473226] ? DS = mfsymbol(mf, D); bestappr(om*op / mfpetersson(DS), 10^8) %3 = 8192/225 ? mf = mfinit([4, 9, -4], 0); ? F = mfeigenbasis(mf)[1]; L = lfunmf(mf, F); ? [v, om] = lfunmfspec(L) %6 = [[1, 10/21, 5/18, 5/24, 5/24, 5/18, 10/21, 1],\ 1.1302643192034974852387822584241400608] ? FS = mfsymbol(mf, F); bestappr(om^2 / mfpetersson(FS), 10^8) %7 = 113246208/325

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 nonnegative 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 parameters [N, k, Vga] of the L-function
defined by `ldata`

, corresponding respectively to
the conductor, the functional equation relating values at s and k-s,
and the gamma shifts of the L-function (see `lfuncreate`

). The gamma
shifts are returned to the current precision.

? L = lfuncreate(1); /* Riemann zeta function */ ? lfunparams(L) %2 = [1, 1, [0]]

The library syntax is `GEN `

.**lfunparams**(GEN ldata, long prec)

Returns the `Ldata`

structure attached to the Θ function
of the lattice attached to the primitive form proportional 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 following computes the Madelung constant:

? L1=lfunqf(matdiagonal([1,1,1])); ? L2=lfunqf(matdiagonal([4,1,1])); ? L3=lfunqf(matdiagonal([4,4,1])); ? F(s)=6*lfun(L2,s)-12*lfun(L3,s)-lfun(L1,s)*(1-8/4^s); ? F(1/2) %5 = -1.7475645946331821906362120355443974035

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)

Creates the Ldata structure (without initialization) corresponding to the
shift of L by d, that is to the function L_{d} such that
L_{d}(s) = L(s-d). If *flag* = 1, return the product L x L_{d} instead.

? Z = lfuncreate(1); \\ zeta(s) ? L = lfunshift(Z,1); \\ zeta(s-1) ? normlp(Vec(lfunlambda(L,s)-lfunlambda(L,3-s))) %3 = 0.E-38 \\ the expansions coincide to 'seriesprecision' ? lfun(L,1) %4 = -0.50000000000000000000000000000000000000 \\ = zeta(0) ? M = lfunshift(Z,1,1); \\ zeta(s)*zeta(s-1) ? normlp(Vec(lfunlambda(M,s)-lfunlambda(M,2-s))) %6 = 2.350988701644575016 E-38 ? lfun(M,2) \\ simple pole at 2, residue zeta(2) %7 = 1.6449340668482264364724151666460251892*x^-1+O(x^0)

The library syntax is `GEN `

.**lfunshift**(GEN L, GEN d, long flag, 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 nonreal 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, long bitprec)

`lim`

being either a positive upper limit or a nonempty real
interval, 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.)

If you pass an `Linit`

to the function, 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)