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 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 expansion 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₁,...,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 expansion 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₀/.), s), for the quadratic Kronecker symbol attached to the fundamental discriminant D₀ = 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₁,...,a_n] such that χ(∏_j g_j^n_j) = 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₁,...,a_n] such that χ(∏_j g_j^n_j) = 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 oo1 oo2
  ? Cf.cyc
  %3 = [6, 2, 2]
  ? Cf.gen
  %4 = [[2, 1; 0, 1], [22, 9; 0, 1], [-6, 7]~]
  ? lfuncreate([Cf, [1,0,0]]); \\  χ(g1) = ζ6, χ(g2) = χ(g3) = 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 CK(𝔪) →  ℝ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 D5
     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 ζ₅ := 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²+(x³+x²+1)y = x²+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τ)^r_m. It is currently assumed that f is a self-dual cuspidal form on Γ₀(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ⁿ), 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₁₀ defined over ℚ(ζ₅), 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₁₀; the third is the factorization of the natural representation of D₁₀ ⊂ S₅; 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₂ |θ(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₁ is the known part of the conductor and D₂ 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(zeta5)
  ? 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(zeta5)
  ? 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 lfuncost0(GEN L, GEN sdom = NULL, long der, long bitprec). Also available is GEN lfuncost(GEN L, GEN dom, long der, long bitprec) when L is 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₁,...,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⁰,..., 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₁,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^k₁+ε), where k₁ = w and even k₁ = 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₁], where k₁ 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₁/x²+r₂/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₁ and L₂ given by L1 and L2. Assume that v_z(L₁) ≥ v_z(L₂) 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}ⁿ η(M_i,1 z)^M_i,2 It is currently assumed that f is a self-dual cuspidal form on Γ₀(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² = F(x) or y²+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).ρ^-Ee^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|^-Ee^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+3N^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₁,r₁],...,[a_n, r_n], [[b₁, R₁],...,[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ⁿ 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).