Since the values of transcendental functions cannot be exactly represented, these functions will always return an inexact object: a real number, a complex number, a p-adic number or a power series. All these objects have a certain finite precision.

As a general rule, which of course in some cases may have exceptions, transcendental functions operate in the following way:

***** If the argument is either a real number or an inexact complex number
(like `1.0 + I`

or `Pi*I`

but not `2 - 3*I`

), then the
computation is done with the precision of the argument.
In the example below, we see that changing the precision to 50 digits does
not matter, because x only had a precision of 19 digits.

? \p 15 realprecision = 19 significant digits (15 digits displayed) ? x = Pi/4 %1 = 0.785398163397448 ? \p 50 realprecision = 57 significant digits (50 digits displayed) ? sin(x) %2 = 0.7071067811865475244

Note that even if the argument is real, the result may be complex (e.g. acos(2.0) or acosh(0.0)). See each individual function help for the definition of the branch cuts and choice of principal value.

***** If the argument is either an integer, a rational, an exact complex
number or a quadratic number, it is first converted to a real
or complex number using the current precision, which can be
view and manipulated using the defaults `realprecision`

(in decimal
digits) or `realbitprecision`

(in bits). This precision can be changed
indifferently

***** in decimal digits: use `\p`

or `default(realprecision,...)`

.

***** in bits: use `\pb`

or `default(realbitprecision,...)`

.

After this conversion, the computation proceeds as above for real or complex arguments.

In library mode, the `realprecision`

does not matter; instead the
precision is taken from the `prec`

parameter which every transcendental
function has. As in `gp`

, this `prec`

is not used when the argument
to a function is already inexact. Note that the argument *prec* stands
for the length in words of a real number, including codewords. Hence we must
have *prec* ≥ 3. (Some functions allow a `bitprec`

argument
instead which allow finer granularity.)

Some accuracies attainable on 32-bit machines cannot be attained
on 64-bit machines for parity reasons. For example the default `gp`

accuracy
is 28 decimal digits on 32-bit machines, corresponding to *prec* having
the value 5, but this cannot be attained on 64-bit machines.

***** If the argument is a polmod (representing an algebraic number),
then the function is evaluated for every possible complex embedding of that
algebraic number. A column vector of results is returned, with one component
for each complex embedding. Therefore, the number of components equals
the degree of the `t_POLMOD`

modulus.

***** If the argument is an intmod or a p-adic, at present only a
few functions like `sqrt`

(square root), `sqr`

(square), `log`

,
`exp`

, powering, `teichmuller`

(Teichmüller character) and
`agm`

(arithmetic-geometric mean) are implemented.

Note that in the case of a 2-adic number, `sqr`

(x) may not be
identical to x*x: for example if x = 1+O(2^5) and y = 1+O(2^5) then
x*y = 1+O(2^5) while `sqr`

(x) = 1+O(2^6). Here, x * x yields the
same result as `sqr`

(x) since the two operands are known to be
*identical*. The same statement holds true for p-adics raised to the
power n, where v_{p}(n) > 0.

**Remark.** If we wanted to be strictly consistent with
the PARI philosophy, we should have x*y = (4 mod 8) and `sqr`

(x) =
(4 mod 32) when both x and y are congruent to 2 modulo 4.
However, since intmod is an exact object, PARI assumes that the modulus
must not change, and the result is hence (0 mod 4) in both cases. On
the other hand, p-adics are not exact objects, hence are treated
differently.

***** If the argument is a polynomial, a power series or a rational function,
it is, if necessary, first converted to a power series using the current
series precision, held in the default `seriesprecision`

. This precision
(the number of significant terms) can be changed using `\ps`

or
`default(seriesprecision,...)`

. Then the Taylor series expansion of the
function around X = 0 (where X is the main variable) is computed to a
number of terms depending on the number of terms of the argument and the
function being computed.

Under `gp`

this again is transparent to the user. When programming in
library mode, however, it is *strongly* advised to perform an explicit
conversion to a power series first, as in

x = gtoser(x, gvar(x), seriesprec)

where the number of significant terms `seriesprec`

can be specified
explicitly. If you do not do this, a global variable `precdl`

is used
instead, to convert polynomials and rational functions to a power series with
a reasonable number of terms; tampering with the value of this global
variable is *deprecated* and strongly discouraged.

***** If the argument is a vector or a matrix, the result is the
componentwise evaluation of the function. In particular, transcendental
functions on square matrices, which are not implemented in the present
version **2.12.0**, will have a different name if they are implemented some day.

The expression x^n is powering.

***** If the exponent n is an integer, then exact operations are performed
using binary (left-shift) powering techniques. By definition, x^0 is
(an empty product interpreted as) an exact 1 in the underlying prime
ring:

? 0.0 ^ 0 %1 = 1 ? (1 + O(2^3)) ^ 0 %2 = 1 ? (1 + O(x)) ^ 0 %3 = 1 ? Mod(2,4)^0 %4 = Mod(1,4) ? Mod(x,x^2)^0 %5 = Mod(1, x^2)

If x is a p-adic number, its precision will increase if v_{p}(n) > 0 and
n != 0. Powering a binary quadratic form (types `t_QFI`

and
`t_QFR`

) returns a representative of the class, which is reduced if the
input was. (In particular, `x^1`

returns x itself, whether it is
reduced or not.)

PARI rewrites the multiplication x * x of two *identical*
objects as x^2. Here, identical means the operands are reference the same
chunk of memory; no equality test is performed. This is no longer true when
more than two arguments are involved.

? a = 1 + O(2); b = a; ? a * a \\ = a^2, precision increases %2 = 1 + O(2^3) ? a * b \\ not rewritten as a^2 %3 = 1 + O(2) ? a*a*a \\ not rewritten as a^3 %4 = 1 + O(2)

***** If the exponent is a rational number p/q, return a solution y of
y^q = x^p if it exists. Beware that this is defined up to q-th roots of 1
in the base field. Intmods modulo composite numbers are not supported.

? 4^(1/2) \\ 4 being a square, this is exact %1 = 2 ? 2^(1/2) \\ now inexact %2 = 1.4142135623730950488016887242096980786 ? Mod(7,19)^(1/2) %3 = Mod(11, 19) \\ is any square root ? sqrt(Mod(7,19)) %4 = Mod(8, 19) \\ is the smallest square root ? Mod(1,4)^(1/2) *** at top-level: Mod(1,4)^(1/2) *** ^ — — *** _^_: not a prime number in gpow: 4.

***** If the exponent is a negative integer or rational number,
an inverse must be computed. For non-invertible `t_INTMOD`

x, this
will fail and (for n an integer) implicitly exhibit a factor of the modulus:

? Mod(4,6)^(-1) *** at top-level: Mod(4,6)^(-1) *** ^ — -- *** _^_: impossible inverse modulo: Mod(2, 6).

Here, a factor 2 is obtained directly. In general, take the gcd of the representative and the modulus. This is most useful when performing complicated operations modulo an integer N whose factorization is unknown. Either the computation succeeds and all is well, or a factor d is discovered and the computation may be restarted modulo d or N/d.

For non-invertible `t_POLMOD`

x, the behavior is the same:

? Mod(x^2, x^3-x)^(-1) *** at top-level: Mod(x^2,x^3-x)^(-1) *** ^ — -- *** _^_: impossible inverse in RgXQ_inv: Mod(x^2, x^3 - x).

Note that the underlying algorihm (subresultant) assumes that the base ring is a domain:

? a = Mod(3*y^3+1, 4); b = y^6+y^5+y^4+y^3+y^2+y+1; c = Mod(a,b); ? c^(-1) *** at top-level: Mod(a,b)^(-1) *** ^ — -- *** _^_: impossible inverse modulo: Mod(2, 4).

In fact c is invertible, but ℤ/4ℤ is not a domain and the algorithm fails. It is possible for the algorithm to succeed in such situations and any returned result will be correct, but chances are that an error will occur first. In this specific case, one should work with 2-adics. In general, one can also try the following approach

? inversemod(a, b) = { my(m, v = variable(b)); m = polsylvestermatrix(polrecip(a), polrecip(b)); m = matinverseimage(m, matid(#m)[,1]); Polrev(m[1..poldegree(b)], v); } ? inversemod(a,b) %2 = Mod(2,4)*y^5 + Mod(3,4)*y^3 + Mod(1,4)*y^2 + Mod(3,4)*y + Mod(2,4)

This is not guaranteed to work either since `matinverseimage`

must also
invert pivots. See Section se:linear_algebra.

For a `t_MAT`

x, the matrix is expected to be square and invertible, except
in the special case `x^(-1)`

which returns a left inverse if one exists
(rectangular x with full column rank).

? x = Mat([1;2]) %1 = [1] [2] ? x^(-1) %2 = [1 0]

***** Finally, if the exponent n is not an rational number, powering is
treated as the transcendental function exp(nlog x), and in particular acts
componentwise on vector or matrices, even square matrices ! (See
Section se:trans.) If x is 0 and n is an inexact 0, this will raise
an exception:

? 4 ^ 1.0 %1 = 4.0000000000000000000000000000000000000 ? 0^ 0.0 *** at top-level: 0^0.0 *** ^ — - *** _^_: domain error in gpow(0,n): n <= 0

The library syntax is `GEN `

for x^n.**gpow**(GEN x, GEN n, long prec)

Catalan's constant G = ∑_{n >= 0}((-1)^n)/((2n+1)^2) = 0.91596....
Note that `Catalan`

is one of the few reserved names which cannot be
used for user variables.

The library syntax is `GEN `

.**mpcatalan**(long prec)

Euler's constant γ = 0.57721.... Note that
`Euler`

is one of the few reserved names which cannot be used for
user variables.

The library syntax is `GEN `

.**mpeuler**(long prec)

The complex number sqrt{-1}.

The library syntax is `GEN `

.**gen _{I}**()

The constant π (3.14159...). Note that `Pi`

is one of the few
reserved names which cannot be used for user variables.

The library syntax is `GEN `

.**mppi**(long prec)

Absolute value of x (modulus if x is complex).
Rational functions are not allowed. Contrary to most transcendental
functions, an exact argument is *not* converted to a real number before
applying `abs`

and an exact result is returned if possible.

? abs(-1) %1 = 1 ? abs(3/7 + 4/7*I) %2 = 5/7 ? abs(1 + I) %3 = 1.414213562373095048801688724

If x is a polynomial, returns -x if the leading coefficient is real and negative else returns x. For a power series, the constant coefficient is considered instead.

The library syntax is `GEN `

.**gabs**(GEN x, long prec)

Principal branch of cos^{-1}(x) = -i log (x + isqrt{1-x^2}).
In particular, Re(acos(x)) ∈ [0,π] and if x ∈ ℝ and |x| > 1,
then acos(x) is complex. The branch cut is in two pieces:
]- oo ,-1] , continuous with quadrant II, and [1,+ oo [, continuous
with quadrant IV. We have acos(x) = π/2 - asin(x) for all
x.

The library syntax is `GEN `

.**gacos**(GEN x, long prec)

Principal branch of cosh^{-1}(x) = 2
log(sqrt{(x+1)/2} + sqrt{(x-1)/2}). In particular,
Re(acosh(x)) ≥ 0 and
Im(acosh(x)) ∈ ]-π,π]; if x ∈ ℝ and x < 1, then
acosh(x) is complex.

The library syntax is `GEN `

.**gacosh**(GEN x, long prec)

Arithmetic-geometric mean of x and y. In the case of complex or negative numbers, the optimal AGM is returned (the largest in absolute value over all choices of the signs of the square roots). p-adic or power series arguments are also allowed. Note that a p-adic agm exists only if x/y is congruent to 1 modulo p (modulo 16 for p = 2). x and y cannot both be vectors or matrices.

The library syntax is `GEN `

.**agm**(GEN x, GEN y, long prec)

Airy [Ai,Bi] functions of argument z.

? [A,B] = airy(1); ? A %2 = 0.13529241631288141552414742351546630617 ? B %3 = 1.2074235949528712594363788170282869954

The library syntax is `GEN `

.**airy**(GEN z, long prec)

Argument of the complex number x, such that -π < arg(x) ≤ π.

The library syntax is `GEN `

.**garg**(GEN x, long prec)

Principal branch of sin^{-1}(x) = -i log(ix + sqrt{1 - x^2}).
In particular, Re(asin(x)) ∈ [-π/2,π/2] and if x ∈ ℝ and
|x| > 1 then asin(x) is complex. The branch cut is in two pieces:
]- oo ,-1], continuous with quadrant II, and [1,+ oo [ continuous
with quadrant IV. The function satisfies i asin(x) =
asinh(ix).

The library syntax is `GEN `

.**gasin**(GEN x, long prec)

Principal branch of sinh^{-1}(x) = log(x + sqrt{1+x^2}). In
particular Im(asinh(x)) ∈ [-π/2,π/2].
The branch cut is in two pieces: ]-i oo ,-i], continuous with quadrant
III and [+i,+i oo [, continuous with quadrant I.

The library syntax is `GEN `

.**gasinh**(GEN x, long prec)

Principal branch of tan^{-1}(x) = log ((1+ix)/(1-ix)) /
2i. In particular the real part of atan(x) belongs to
]-π/2,π/2[.
The branch cut is in two pieces:
]-i oo ,-i[, continuous with quadrant IV, and ]i,+i oo [ continuous
with quadrant II. The function satisfies atan(x) =
-iatanh(ix) for all x != ± i.

The library syntax is `GEN `

.**gatan**(GEN x, long prec)

Principal branch of tanh^{-1}(x) = log ((1+x)/(1-x)) / 2. In
particular the imaginary part of atanh(x) belongs to
[-π/2,π/2]; if x ∈ ℝ and |x| > 1 then atanh(x) is complex.

The library syntax is `GEN `

.**gatanh**(GEN x, long prec)

Bernoulli number B_{n},
where B_{0} = 1, B_{1} = -1/2, B_{2} = 1/6,..., expressed as a rational number.
The argument n should be a non-negative integer. The function `bervec`

creates a cache of successive Bernoulli numbers which greatly speeds up
later calls to `bernfrac`

:

? bernfrac(20000); time = 107 ms. ? bernvec(10000); \\ cache B_{0}, B_{2}, ..., B_20000 time = 35,957 ms. ? bernfrac(20000); \\ now instantaneous ?

The library syntax is `GEN `

.**bernfrac**(long n)

Bernoulli polynomial B_{n} in variable v.

? bernpol(1) %1 = x - 1/2 ? bernpol(3) %2 = x^3 - 3/2*x^2 + 1/2*x

The library syntax is `GEN `

where **bernpol**(long n, long v = -1)`v`

is a variable number.

Bernoulli number
B_{n}, as `bernfrac`

, but B_{n} is returned as a real number
(with the current precision). The argument n should be a non-negative
integer. The function slows down as the precision increases:

? \p1000 ? bernreal(200000); time = 5 ms. ? \p10000 ? bernreal(200000); time = 18 ms. ? \p100000 ? bernreal(200000); time = 84 ms.

The library syntax is `GEN `

.**bernreal**(long n, long prec)

Returns a vector containing, as rational numbers,
the Bernoulli numbers B_{0}, B_{2},..., B_{2n}:

? bernvec(5) \\ B_{0}, B_{2}..., B_10 %1 = [1, 1/6, -1/30, 1/42, -1/30, 5/66] ? bernfrac(10) %2 = 5/66

This routine uses a lot of memory but is much faster than
repeated calls to `bernfrac`

:

? forstep(n = 2, 10000, 2, bernfrac(n)) time = 41,522 ms. ? bernvec(5000); time = 4,784 ms.

The computed Bernoulli numbers are stored in an incremental
cache which makes later calls to `bernfrac`

and `bernreal`

instantaneous in the cache range: re-running the same previous `bernfrac`

s
after the `bernvec`

call gives:

? forstep(n = 2, 10000, 2, bernfrac(n)) time = 1 ms.

The time and space complexity of this function are
Õ(n^2); in the feasible range n ≤ 10^5 (requires about 2 hours),
the practical time complexity is closer to Õ(n^{log2 6}).

The library syntax is `GEN `

.**bernvec**(long n)

H^1-Bessel function of index *nu* and argument x.

The library syntax is `GEN `

.**hbessel1**(GEN nu, GEN x, long prec)

H^2-Bessel function of index *nu* and argument x.

The library syntax is `GEN `

.**hbessel2**(GEN nu, GEN x, long prec)

I-Bessel function of index *nu* and
argument x. If x converts to a power series, the initial factor
(x/2)^ν/Γ(ν+1) is omitted (since it cannot be represented in PARI
when ν is not integral).

The library syntax is `GEN `

.**ibessel**(GEN nu, GEN x, long prec)

J-Bessel function of index *nu* and
argument x. If x converts to a power series, the initial factor
(x/2)^ν/Γ(ν+1) is omitted (since it cannot be represented in PARI
when ν is not integral).

The library syntax is `GEN `

.**jbessel**(GEN nu, GEN x, long prec)

J-Bessel function of half integral index.
More precisely, `besseljh`

(n,x) computes J_{n+1/2}(x) where n
must be of type integer, and x is any element of ℂ. In the
present version **2.12.0**, this function is not very accurate when x is small.

The library syntax is `GEN `

.**jbesselh**(GEN n, GEN x, long prec)

K-Bessel function of index *nu* and argument x.

The library syntax is `GEN `

.**kbessel**(GEN nu, GEN x, long prec)

Deprecated alias for `bessely`

.

The library syntax is `GEN `

.**nbessel**(GEN nu, GEN x, long prec)

Y-Bessel function of index *nu* and argument x.

The library syntax is `GEN `

.**nbessel**(GEN nu, GEN x, long prec)

Cosine of x.

The library syntax is `GEN `

.**gcos**(GEN x, long prec)

Hyperbolic cosine of x.

The library syntax is `GEN `

.**gcosh**(GEN x, long prec)

Cotangent of x.

The library syntax is `GEN `

.**gcotan**(GEN x, long prec)

Hyperbolic cotangent of x.

The library syntax is `GEN `

.**gcotanh**(GEN x, long prec)

Principal branch of the dilogarithm of x,
i.e. analytic continuation of the power series log_{2}(x) = ∑_{n ≥ 1}x^n/n^2.

The library syntax is `GEN `

.**dilog**(GEN x, long prec)

Exponential integral ∫_{x}^ oo (e^{-t})/(t)dt =
`incgam`

(0, x), where the latter expression extends the function
definition from real x > 0 to all complex x != 0.

If n is present, we must have x > 0; the function returns the
n-dimensional vector [`eint1`

(x),...,`eint1`

(nx)]. Contrary to
other transcendental functions, and to the default case (n omitted), the
values are correct up to a bounded *absolute*, rather than relative,
error 10^{-n}, where n is `precision`

(x) if x is a `t_REAL`

and defaults to `realprecision`

otherwise. (In the most important
application, to the computation of L-functions via approximate functional
equations, those values appear as weights in long sums and small individual
relative errors are less useful than controlling the absolute error.) This is
faster than repeatedly calling `eint1(i * x)`

, but less precise.

The library syntax is `GEN `

.
Also available is **veceint1**(GEN x, GEN n = NULL, long prec)`GEN `

.**eint1**(GEN x, long prec)

Complementary error function, analytic continuation of
(2/sqrtπ)∫_{x}^ oo e^{-t^2}dt = `incgam`

(1/2,x^2)/sqrtπ,
where the latter expression extends the function definition from real x to
all complex x != 0.

The library syntax is `GEN `

.**gerfc**(GEN x, long prec)

Variants of Dedekind's η function.
If *flag* = 0, return ∏_{n = 1}^ oo (1-q^n), where q depends on x
in the following way:

***** q = e^{2iπ x} if x is a *complex number* (which must then
have positive imaginary part); notice that the factor q^{1/24} is
missing!

***** q = x if x is a `t_PADIC`

, or can be converted to a
*power series* (which must then have positive valuation).

If *flag* is non-zero, x is converted to a complex number and we return the
true η function, q^{1/24}∏_{n = 1}^ oo (1-q^n),
where q = e^{2iπ x}.

The library syntax is `GEN `

.**eta0**(GEN z, long flag, long prec)

Also available is `GEN `

(**trueeta**(GEN x, long prec)*flag* = 1).

Exponential of x. p-adic arguments with positive valuation are accepted.

The library syntax is `GEN `

.
For a **gexp**(GEN x, long prec)`t_PADIC`

x, the function
`GEN `

is also available.**Qp_exp**(GEN x)

Return exp(x)-1, computed in a way that is also accurate
when the real part of x is near 0.
A naive direct computation would suffer from catastrophic cancellation;
PARI's direct computation of exp(x) alleviates this well known problem at
the expense of computing exp(x) to a higher accuracy when x is small.
Using `expm1`

is recommended instead:

? default(realprecision, 10000); x = 1e-100; ? a = expm1(x); time = 4 ms. ? b = exp(x)-1; time = 4 ms. ? default(realprecision, 10040); x = 1e-100; ? c = expm1(x); \\ reference point ? abs(a-c)/c \\ relative error in expm1(x) %7 = 1.4027986153764843997 E-10019 ? abs(b-c)/c \\ relative error in exp(x)-1 %8 = 1.7907031188259675794 E-9919

As the example above shows, when x is near 0,
`expm1`

is more accurate than `exp(x)-1`

.

The library syntax is `GEN `

.**gexpm1**(GEN x, long prec)

For s a complex number, evaluates Euler's gamma
function
Γ(s) = ∫_{0}^ oo t^{s-1}exp(-t)dt.
Error if s is a non-positive integer, where Γ has a pole.

For s a `t_PADIC`

, evaluates the Morita gamma function at s, that
is the unique continuous p-adic function on the p-adic integers
extending Γ_{p}(k) = (-1)^k ∏_{j < k}'j, where the prime means that p
does not divide j.

? gamma(1/4 + O(5^10)) %1= 1 + 4*5 + 3*5^4 + 5^6 + 5^7 + 4*5^9 + O(5^10) ? algdep(%,4) %2 = x^4 + 4*x^2 + 5

The library syntax is `GEN `

.
For a **ggamma**(GEN s, long prec)`t_PADIC`

x, the function `GEN `

is
also available.**Qp_gamma**(GEN x)

Gamma function evaluated at the argument x+1/2.

The library syntax is `GEN `

.**ggammah**(GEN x, long prec)

Returns the value at t of the inverse Mellin transform
G initialized by `gammamellininvinit`

.

? G = gammamellininvinit([0]); ? gammamellininv(G, 2) - 2*exp(-Pi*2^2) %2 = -4.484155085839414627 E-44

The alternative shortcut

gammamellininv(A,t,m)

for

gammamellininv(gammamellininvinit(A,m), t)

is available.

The library syntax is `GEN `

.**gammamellininv**(GEN G, GEN t, long m, long bitprec)

Return the first n terms of the asymptotic expansion at infinity
of the m-th derivative K^{(m)}(t) of the inverse Mellin transform of the
function
f(s) = Γ_ℝ(s+a_{1})...Γ_ℝ(s+a_{d}) ,
where `A`

is the vector [a_{1},...,a_{d}] and
Γ_ℝ(s) = π^{-s/2} Γ(s/2) (Euler's `gamma`

).
The result is a vector
[M[1]...M[n]] with M[1] = 1, such that
K^{(m)}(t) = sqrt{2^{d+1}/d}t^{a+m(2/d-1)}e^{-dπ t^{2/d}}
∑_{n ≥ 0} M[n+1] (π t^{2/d})^{-n}
with a = (1-d+∑_{1 ≤ j ≤ d}a_{j})/d.

The library syntax is `GEN `

.**gammamellininvasymp**(GEN A, long precdl, long n)

Initialize data for the computation by `gammamellininv`

of
the m-th derivative of the inverse Mellin transform of the function
f(s) = Γ_ℝ(s+a_{1})...Γ_ℝ(s+a_{d})
where `A`

is the vector [a_{1},...,a_{d}] and
Γ_ℝ(s) = π^{-s/2} Γ(s/2) (Euler's `gamma`

). This is the
special case of Meijer's G functions used to compute L-values via the
approximate functional equation.

**Caveat.** Contrary to the PARI convention, this function
guarantees an *absolute* (rather than relative) error bound.

For instance, the inverse Mellin transform of Γ_ℝ(s) is 2exp(-π z^2):

? G = gammamellininvinit([0]); ? gammamellininv(G, 2) - 2*exp(-Pi*2^2) %2 = -4.484155085839414627 E-44

The inverse Mellin transform of Γ_ℝ(s+1) is 2 zexp(-π z^2), and its second derivative is 4π z exp(-π z^2)(2π z^2 - 3):

? G = gammamellininvinit([1], 2); ? a(z) = 4*Pi*z*exp(-Pi*z^2)*(2*Pi*z^2-3); ? b(z) = gammamellininv(G,z); ? t(z) = b(z) - a(z); ? t(3/2) %3 = -1.4693679385278593850 E-39

The library syntax is `GEN `

.**gammamellininvinit**(GEN A, long m, long bitprec)

General hypergeometric function, where `N`

and `D`

are
the vector of parameters in the numerator and denominator respectively,
evaluated at the complex argument z.

This function implements hypergeometric functions
_pF_{q}((a_{i})_{1 ≤ i ≤ p},(b_{j})_{1 ≤ j ≤ q};z)
= ∑_{n ≥ 0}(∏_{1 ≤ i ≤ p}(a_{i})_{n})/(∏_{1 ≤ j ≤ q}(b_{j})_{n})
(z^n)/(n!) ,
where (a)_{n} = a(a+1)...(a+n-1) is the rising Pochammer symbol. For this
to make sense, none of the b_{j} must be a negative or zero integer. The
corresponding general GP command is

hypergeom([a1,a2,...,ap], [b1,b2,...,bq], z)

Whenever p = 1 or q = 1, a one-element vector can be
replaced by the element it contains. Whenever p = 0 or q = 0, an empty
vector can be omitted. For instance hypergeom(,b,z) computes _0F_{1}(;b;z).

We distinguish three kinds of such functions according to their radius of convergence R:

***** q ≥ p: R = oo .

***** q = p-1: R = 1. Nonetheless, by integral representations, _pF_{q}
can be analytically continued outside the disc of convergence.

***** q ≤ p-2: R = 0. By integral representations, one can make sense of
the function in a suitable domain.

The list of implemented functions and their domain of validity in our implementation is as follows:

`F01`

: `hypergeom(,a,z)`

(or `[a]`

).
This is essentially a Bessel function and computed as such. R = oo .

`F21`

: `hypergeom([a,b],c,z)`

(or `[c]`

).
R = 1, extended by
(Γ(c))/(Γ(b)Γ(c-b))
∫_{0}^1 t^{b-1}(1-t)^{c-b-1}(1-zt)^adt .
This is Gauss's Hypergeometric function, and almost all of the implementation
work is done for this function.

`F32`

: `hypergeom([a,b,c],[d,e],z)`

. R = 1, extended by
(Γ(e))/(Γ(c)Γ(e-c))
∫_{0}^1t^{c-1}(1-t)^{e-c-1}_2F_{1}(a,b;d;tz)dt .

`F20`

: `hypergeom([a,b],,z)`

. R = 0, computed as
(1)/(Γ(a))∫_{0}^ oo t^{a-1}(1-zt)^{-b}e^{-t}dt .

`F31`

: `hypergeom([a,b,c],d,z)`

(or `[d]`

). R = 0, computed as
(1)/(Γ(a))∫_{0}^ oo t^{a-1}e^{-t}_2F_{1}(b,c;d;tz)dt .

For other inputs: if R = oo or R = 1 and |z| < 1- ϵ is not too close to the circle of convergence, we simply sum the series.

? hypergeom([3,2], 3.4, 0.7) \\ 2F1(3,2; 3.4; 0.7) %1 = 7.9999999999999999999999999999999999999 ? a=5/3; T1=hypergeom([1,1,1],[a,a],1) \\ 3F2(1,1,1; a,a; 1) %2 = 3.1958592952314032651578713968927593818 ? T2=hypergeom([2,1,1],[a+1,a+1],1) %3 = 1.6752931349345765309211012564734179541 ? T3=hypergeom([2*a-1,1,1],[a+1,a+1],1) %4 = 1.9721037126267142061807688820853354440 ? T1 + (a-1)^2/(a^2*(2*a-3)) * (T2-2*(a-1)*T3) \\ - gamma(a)^2/((2*a-3)*gamma(2*a-2)) %5 = -1.880790961315660013 E-37 \\ ~ 0

This identity is due to Bercu.

The library syntax is `GEN `

.**hypergeom**(GEN N = NULL, GEN D = NULL, GEN z, long prec)

U-confluent hypergeometric function with complex
parameters a, b, z. Note that _2F_{0}(a,b,z) = (-z)^{-a}U(a, a+1-b, -1/z),

? hyperu(1, 3/2, I) %1 = 0.23219... - 0.80952...*I ? -I * hypergeom([1, 1+1-3/2], [], -1/I) %2 = 0.23219... - 0.80952...*I

The library syntax is `GEN `

.**hyperu**(GEN a, GEN b, GEN z, long prec)

Incomplete gamma function ∫_{x}^ oo e^{-t}t^{s-1}dt, extended by
analytic continuation to all complex x, s not both 0. The relative error
is bounded in terms of the precision of s (the accuracy of x is ignored
when determining the output precision). When g is given, assume that
g = Γ(s). For small |x|, this will speed up the computation.

The library syntax is `GEN `

.
Also available is **incgam0**(GEN s, GEN x, GEN g = NULL, long prec)`GEN `

.**incgam**(GEN s, GEN x, long prec)

Complementary incomplete gamma function.
The arguments x and s are complex numbers such that s is not a pole of
Γ and |x|/(|s|+1) is not much larger than 1 (otherwise the
convergence is very slow). The result returned is ∫_{0}^x
e^{-t}t^{s-1}dt.

The library syntax is `GEN `

.**incgamc**(GEN s, GEN x, long prec)

Lambert W function, solution of the implicit equation xe^x = y, for y > 0.

The library syntax is `GEN `

.**glambertW**(GEN y, long prec)

Principal branch of the logarithm of the gamma function of x. This
function is analytic on the complex plane with non-positive integers
removed, and can have much larger arguments than `gamma`

itself.

For x a power series such that x(0) is not a pole of `gamma`

,
compute the Taylor expansion. (PARI only knows about regular power series
and can't include logarithmic terms.)

? lngamma(1+x+O(x^2)) %1 = -0.57721566490153286060651209008240243104*x + O(x^2) ? lngamma(x+O(x^2)) *** at top-level: lngamma(x+O(x^2)) *** ^ — — — — — -- *** lngamma: domain error in lngamma: valuation != 0 ? lngamma(-1+x+O(x^2)) *** lngamma: Warning: normalizing a series with 0 leading term. *** at top-level: lngamma(-1+x+O(x^2)) *** ^ — — — — — — -- *** lngamma: domain error in intformal: residue(series, pole) != 0

The library syntax is `GEN `

.**glngamma**(GEN x, long prec)

Principal branch of the natural logarithm of
x ∈ ℂ^{*}, i.e. such that Im(log(x)) ∈ ]-π,π].
The branch cut lies
along the negative real axis, continuous with quadrant 2, i.e. such that
lim_{b → 0^+} log (a+bi) = log a for a ∈ ℝ^{*}. The result is complex
(with imaginary part equal to π) if x ∈ ℝ and x < 0. In general,
the algorithm uses the formula
log(x) ~ (π)/(2agm(1, 4/s)) - m log 2,
if s = x 2^m is large enough. (The result is exact to B bits provided
s > 2^{B/2}.) At low accuracies, the series expansion near 1 is used.

p-adic arguments are also accepted for x, with the convention that log(p) = 0. Hence in particular exp(log(x))/x is not in general equal to 1 but to a (p-1)-th root of unity (or ±1 if p = 2) times a power of p.

The library syntax is `GEN `

.
For a **glog**(GEN x, long prec)`t_PADIC`

x, the function
`GEN `

is also available.**Qp_log**(GEN x)

Return log(1+x), computed in a way that is also accurate
when the real part of x is near 0. This is the reciprocal function
of `expm1`

(x) = exp(x)-1.

? default(realprecision, 10000); x = Pi*1e-100; ? (expm1(log1p(x)) - x) / x %2 = -7.668242895059371866 E-10019 ? (log1p(expm1(x)) - x) / x %3 = -7.668242895059371866 E-10019

When x is small, this function is both faster and more accurate than log(1+x):

? \p38 ? x = 1e-20; ? localprec(100); c = log1p(x); \\ reference point ? a = log1p(x); abs((a - c)/c) %6 = 0.E-38 ? b = log(1+x); abs((b - c)/c) \\ slightly less accurate %7 = 1.5930919111324522770 E-38 ? for (i=1,10^5,log1p(x)) time = 81 ms. ? for (i=1,10^5,log(1+x)) time = 100 ms. \\ slower, too

The library syntax is `GEN `

.**glog1p**(GEN x, long prec)

One of the different polylogarithms, depending on *flag*:

If *flag* = 0 or is omitted: m-th polylogarithm of x, i.e. analytic
continuation of the power series Li_{m}(x) = ∑_{n ≥ 1}x^n/n^m
(x < 1). Uses the functional equation linking the values at x and 1/x
to restrict to the case |x| ≤ 1, then the power series when
|x|^2 ≤ 1/2, and the power series expansion in log(x) otherwise.

Using *flag*, computes a modified m-th polylogarithm of x.
We use Zagier's notations; let Re_{m} denote Re or Im depending
on whether m is odd or even:

If *flag* = 1: compute ~ D_{m}(x), defined for |x| ≤ 1 by
Re_{m}(∑_{k = 0}^{m-1} ((-log|x|)^k)/(k!)Li_{m-k}(x)
+((-log|x|)^{m-1})/(m!)log|1-x|).

If *flag* = 2: compute D_{m}(x), defined for |x| ≤ 1 by
Re_{m}(∑_{k = 0}^{m-1}((-log|x|)^k)/(k!)Li_{m-k}(x)
-(1)/(2)((-log|x|)^m)/(m!)).

If *flag* = 3: compute P_{m}(x), defined for |x| ≤ 1 by
Re_{m}(∑_{k = 0}^{m-1}(2^kB_{k})/(k!)(log|x|)^kLi_{m-k}(x)
-(2^{m-1}B_{m})/(m!)(log|x|)^m).

These three functions satisfy the functional equation
f_{m}(1/x) = (-1)^{m-1}f_{m}(x).

The library syntax is `GEN `

.
Also available is
**polylog0**(long m, GEN x, long flag, long prec)`GEN `

(**gpolylog**(long m, GEN x, long prec)*flag* = 0).

The ψ-function of x, i.e. the logarithmic derivative Γ'(x)/Γ(x).

The library syntax is `GEN `

.**gpsi**(GEN x, long prec)

Sine of x.

The library syntax is `GEN `

.**gsin**(GEN x, long prec)

Cardinal sine of x, i.e. sin(x)/x if x != 0, 1 otherwise.
Note that this function also allows to compute
(1-cos(x)) / x^2 = `sinc`

(x/2)^2 / 2
accurately near x = 0.

The library syntax is `GEN `

.**gsinc**(GEN x, long prec)

Hyperbolic sine of x.

The library syntax is `GEN `

.**gsinh**(GEN x, long prec)

Square of x. This operation is not completely straightforward, i.e. identical to x * x, since it can usually be computed more efficiently (roughly one-half of the elementary multiplications can be saved). Also, squaring a 2-adic number increases its precision. For example,

? (1 + O(2^4))^2 %1 = 1 + O(2^5) ? (1 + O(2^4)) * (1 + O(2^4)) %2 = 1 + O(2^4)

Note that this function is also called whenever one multiplies two objects
which are known to be *identical*, e.g. they are the value of the same
variable, or we are computing a power.

? x = (1 + O(2^4)); x * x %3 = 1 + O(2^5) ? (1 + O(2^4))^4 %4 = 1 + O(2^6)

(note the difference between `%2`

and `%3`

above).

The library syntax is `GEN `

.**gsqr**(GEN x)

Principal branch of the square root of x, defined as sqrt{x} = exp(log x / 2). In particular, we have Arg(sqrt(x)) ∈ ]-π/2, π/2], and if x ∈ ℝ and x < 0, then the result is complex with positive imaginary part.

Intmod a prime p, `t_PADIC`

and `t_FFELT`

are allowed as arguments. In
the first 2 cases (`t_INTMOD`

, `t_PADIC`

), the square root (if it
exists) which is returned is the one whose first p-adic digit is in the
interval [0,p/2]. For other arguments, the result is undefined.

The library syntax is `GEN `

.
For a **gsqrt**(GEN x, long prec)`t_PADIC`

x, the function
`GEN `

is also available.**Qp_sqrt**(GEN x)

Principal branch of the nth root of x, i.e. such that Arg(sqrtn(x)) ∈ ]-π/n, π/n]. Intmod a prime and p-adics are allowed as arguments.

If z is present, it is set to a suitable root of unity allowing to recover all the other roots. If it was not possible, z is set to zero. In the case this argument is present and no nth root exist, 0 is returned instead of raising an error.

? sqrtn(Mod(2,7), 2) %1 = Mod(3, 7) ? sqrtn(Mod(2,7), 2, &z); z %2 = Mod(6, 7) ? sqrtn(Mod(2,7), 3) *** at top-level: sqrtn(Mod(2,7),3) *** ^ — — — — — -- *** sqrtn: nth-root does not exist in gsqrtn. ? sqrtn(Mod(2,7), 3, &z) %2 = 0 ? z %3 = 0

The following script computes all roots in all possible cases:

sqrtnall(x,n)= { my(V,r,z,r2); r = sqrtn(x,n, &z); if (!z, error("Impossible case in sqrtn")); if (type(x) == "t_INTMOD" || type(x)=="t_PADIC", r2 = r*z; n = 1; while (r2!=r, r2*=z;n++)); V = vector(n); V[1] = r; for(i=2, n, V[i] = V[i-1]*z); V } addhelp(sqrtnall,"sqrtnall(x,n):compute the vector of nth-roots of x");

The library syntax is `GEN `

.
If x is a **gsqrtn**(GEN x, GEN n, GEN *z = NULL, long prec)`t_PADIC`

, the function
`GEN `

is also available.**Qp_sqrtn**(GEN x, GEN n, GEN *z)

Tangent of x.

The library syntax is `GEN `

.**gtan**(GEN x, long prec)

Hyperbolic tangent of x.

The library syntax is `GEN `

.**gtanh**(GEN x, long prec)

Teichmüller character of the p-adic number x, i.e. the unique
(p-1)-th root of unity congruent to x / p^{vp(x)} modulo p.
If x is of the form [p,n], for a prime p and integer n,
return the lifts to ℤ of the images of i + O(p^n) for
i = 1,..., p-1, i.e. all roots of 1 ordered by residue class modulo
p. Such a vector can be fed back to `teichmuller`

, as the
optional argument `tab`

, to speed up later computations.

? z = teichmuller(2 + O(101^5)) %1 = 2 + 83*101 + 18*101^2 + 69*101^3 + 62*101^4 + O(101^5) ? z^100 %2 = 1 + O(101^5) ? T = teichmuller([101, 5]); ? teichmuller(2 + O(101^5), T) %4 = 2 + 83*101 + 18*101^2 + 69*101^3 + 62*101^4 + O(101^5)

As a rule of thumb, if more than
p / 2(log_{2}(p) + `hammingweight`

(p))
values of `teichmuller`

are to be computed, then it is worthwile to
initialize:

? p = 101; n = 100; T = teichmuller([p,n]); \\ instantaneous ? for(i=1,10^3, vector(p-1, i, teichmuller(i+O(p^n), T))) time = 60 ms. ? for(i=1,10^3, vector(p-1, i, teichmuller(i+O(p^n)))) time = 1,293 ms. ? 1 + 2*(log(p)/log(2) + hammingweight(p)) %8 = 22.316[...]

Here the precompuation induces a speedup by a factor 1293/ 60 ~ 21.5.

**Caveat.**
If the accuracy of `tab`

(the argument n above) is lower than the
precision of x, the *former* is used, i.e. the cached value is not
refined to higher accuracy. It the accuracy of `tab`

is larger, then
the precision of x is used:

? Tlow = teichmuller([101, 2]); \\ lower accuracy ! ? teichmuller(2 + O(101^5), Tlow) %10 = 2 + 83*101 + O(101^5) \\ no longer a root of 1 ? Thigh = teichmuller([101, 10]); \\ higher accuracy ? teichmuller(2 + O(101^5), Thigh) %12 = 2 + 83*101 + 18*101^2 + 69*101^3 + 62*101^4 + O(101^5)

The library syntax is `GEN `

.**teichmuller**(GEN x, GEN tab = NULL)

Also available are the functions `GEN `

(**teich**(GEN x)`tab`

is
`NULL`

) as well as
`GEN `

.**teichmullerinit**(long p, long n)

Jacobi sine theta-function
θ_{1}(z, q) = 2q^{1/4} ∑_{n ≥ 0} (-1)^n q^{n(n+1)} sin((2n+1)z).

The library syntax is `GEN `

.**theta**(GEN q, GEN z, long prec)

k-th derivative at z = 0 of `theta`

(q,z).

The library syntax is `GEN `

.**thetanullk**(GEN q, long k, long prec)

`GEN `

returns the vector
of all (d^iθ)/(dz^i)(q,0) for all odd i = 1, 3,..., 2k-1.
**vecthetanullk**(GEN q, long k, long prec)`GEN `

returns
**vecthetanullk_tau**(GEN tau, long k, long prec)`vecthetanullk_tau`

at q = exp(2iπ `tau`

).

One of Weber's three f functions.
If *flag* = 0, returns
f(x) = exp(-iπ/24).η((x+1)/2)/η(x) {such that}
j = (f^{24}-16)^3/f^{24},
where j is the elliptic j-invariant (see the function `ellj`

).
If *flag* = 1, returns
f_{1}(x) = η(x/2)/η(x) {such that}
j = (f_{1}^{24}+16)^3/f_{1}^{24}.
Finally, if *flag* = 2, returns
f_{2}(x) = sqrt{2}η(2x)/η(x) {such that}
j = (f_{2}^{24}+16)^3/f_{2}^{24}.
Note the identities f^8 = f_{1}^8+f_{2}^8 and ff_1f_{2} = sqrt2.

The library syntax is `GEN `

.
Also available are **weber0**(GEN x, long flag, long prec)`GEN `

,
**weberf**(GEN x, long prec)`GEN `

and **weberf1**(GEN x, long prec)`GEN `

.**weberf2**(GEN x, long prec)

For s != 1 a complex number, Riemann's zeta
function ζ(s) = ∑_{n ≥ 1}n^{-s},
computed using the Euler-Maclaurin summation formula, except
when s is of type integer, in which case it is computed using
Bernoulli numbers for s ≤ 0 or s > 0 and
even, and using modular forms for s > 0 and odd. Power series
are also allowed:

? zeta(2) - Pi^2/6 %1 = 0.E-38 ? zeta(1+x+O(x^3)) %2 = 1.0000000000000000000000000000000000000*x^-1 + \ 0.57721566490153286060651209008240243104 + O(x)

For s != 1 a p-adic number, Kubota-Leopoldt zeta function at s, that
is the unique continuous p-adic function on the p-adic integers
that interpolates the values of (1 - p^{-k}) ζ(k) at negative
integers k such that k = 1 (mod p-1) (resp. k is odd) if
p is odd (resp. p = 2). Power series are not allowed in this case.

? zeta(-3+O(5^10)) %1 = 4*5^-1 + 4 + 3*5 + 4*5^3 + 4*5^5 + 4*5^7 + O(5^9))))) ? (1-5^3) * zeta(-3) %2 = -1.0333333333333333333333333333333333333 ? bestappr(%) %3 = -31/30 ? zeta(-3+O(5^10)) - (-31/30) %4 = O(5^9)

The library syntax is `GEN `

.**gzeta**(GEN s, long prec)

Hurwitz zeta function ζ(s,x) = ∑_{n ≥ 0}(n+x)^{-s} and
analytically continued, with s != 1 and x not a negative or zero
integer. Note that ζ(s,1) = ζ(s). s can also be a polynomial,
rational function, or power series. If `der`

is positive, compute the
`der`

'th derivative with respect to s. Note that the derivative
with respect to x is simply -sζ(s+1,x).

? zetahurwitz(Pi,Pi) %1 = 0.056155444497585099925180502385781494484 ? zetahurwitz(2,1) - zeta(2) %2 = -2.350988701644575016 E-38 ? zetahurwitz(Pi,3) - (zeta(Pi)-1-1/2^Pi) %3 = -2.2040519077917890774 E-39 ? zetahurwitz(-7/2,1) - zeta(-7/2) %4 = -2.295887403949780289 E-41 ? zetahurwitz(-2.3,Pi+I*log(2)) %5 = -5.1928369229555125820137832704455696057\ - 6.1349660138824147237884128986232049582*I ? zetahurwitz(-1+x^2+O(x^3),1) %6 = -0.083333333333333333333333333333333333333\ - 0.16542114370045092921391966024278064276*x^2 + O(x^3) ? zetahurwitz(1+x+O(x^4),2) %7 = 1.0000000000000000000000000000000000000*x^-1\ - 0.42278433509846713939348790991759756896\ + 0.072815845483676724860586375874901319138*x + O(x^2) ? zetahurwitz(2,1,2) \\ zeta''(2) %8 = 1.9892802342989010234208586874215163815

The library syntax is `GEN `

.**zetahurwitz**(GEN s, GEN x, long der, long bitprec)

For s a vector of positive integers such that s[1] ≥ 2,
returns the multiple zeta value (MZV)
ζ(s_{1},..., s_{k}) = ∑_{n1 > ... > nk > 0} n_{1}^{-s1}...n_{k}^{-sk}.

? zetamult([2,1]) - zeta(3) \\ Euler's identity %1 = 0.E-38

If the bit precision is B, this function runs in time Õ(k B^2).

If T is provided, it must be the output of `zetamultinit`

(w) for
some w ≥ s_{1} +...+ s_{k} and will provide a small speed up,
usually about 10%.

? T = zetamultinit(20); s = [2,1,1,1,1,1,1,1,1]; ? for(i=1,10^3, zetamult(s)) time = 373 ms. ? for(i=1,10^3, zetamult(s, T)) \\ faster time = 279 ms. ? zetamult(vector(10,i,2), T) %4 = 1.7165384749821433018378232207719985786 E-10 ? zetamult(vector(11,i,2), T) \\ overshoot *** at top-level: zetamult(vector(11,i *** ^ — — — — — — -- *** zetamult: domain error in zetamult: weight > 20

In addition to the above format (`avec`

), the function
also accepts an internal binary format `evec`

(each s_{i} is replaced
by s_{i} bits, all of them 0 but the last one), and an `index`

format
(if e is the positive integer attached the `evec`

vector of
bits, the index is the integer e + 2^{k-2}). The function
`zetamultconvert`

allows to pass from one format to the other; the
function `zetamultall`

computes simultaneously all MZVs of weight
∑_{i ≤ k} s_{i} up to n.

The library syntax is `GEN `

.
Also available is **zetamult0**(GEN s, GEN T = NULL, long prec)`GEN `

.**zetamult**(GEN avec, long prec)

List of all multiple zeta values for weight s_{1} +...+ s_{k} up to n.
The function returns a vector with 2^{n-1}-1 components whose i-th entry
is the MZV of `index`

i (see `zetamult`

).

? z = zetamultall(5); ? z[10] %2 = 0.22881039760335375976874614894168879193 ? zetamultconvert(10) \\ convert index 10 to avec %3 = Vecsmall([3, 2]) ? zetamult(%) %4 = 0.22881039760335375976874614894168879193 ? zetamult(10) %5 = 0.22881039760335375976874614894168879193

If the bit precision is B, this function runs in time O(2^n n B^2) for an output of size O(2^n B).

The library syntax is `GEN `

.**zetamultall**(long n, long prec)

`a`

being either an `evec`

, `avec`

, or index `m`

,
converts into `evec`

(`fl = 0`

), `avec`

(`fl = 1`

, default), or
index `m`

(`fl = 2`

); see `zetamult`

for explanations.

? zetamultconvert(10) %1 = Vecsmall([3, 2]) ? zetamultconvert(13) %2 = Vecsmall([2, 2, 1]) ? zetamultconvert(10, 0) %3 = Vecsmall([0, 0, 1, 0, 1]) ? zetamultconvert(13, 0) %4 = Vecsmall([0, 1, 0, 1, 1])

The last two lines imply that [3,2] and [2,2,1]
are dual (reverse order of bits and swap 0 and 1 in `evec`

form).
Hence they have the same zeta value:

? zetamult([3,2]) %5 = 0.22881039760335375976874614894168879193 ? zetamult([2,2,1]) %6 = 0.22881039760335375976874614894168879193

The library syntax is `GEN `

.**zetamultconvert**(GEN a, long fl)

Initialize data (depending on the precision) used to compute
multiple zeta values at
integral points s = [s_{1},...,s_{k}] for any s_{1} +...+ s_{k} ≤
*maxw*. The corresponding data is inexpensive to compute or store
and provides a small speedup (usually about 10%) when multiple values are to
be computed at a given accuracy.

? for(i = 1, 2^12-1, zetamult(i)) time = 1,413 ms ? T = zetamultinit(13); \\ instantaneous ? for(i = 1, 2^12-1, zetamult(i, T)) \\ used cached data time = 1,315 ms ? zetamultall(12); \\ much faster ! time = 27 ms ? T=zetamultinit(102); sizebyte(T) \\ small even for huge weights time = 5 ms. %5 = 1440504 ? for(i = 1, 2^5, zetamult(2^100+i)) time = 633 ms. ? for(i = 1, 2^5, zetamult(2^100+i, T)) time = 550 ms.

For small weights, `zetamultall`

will be much more efficient; but it
is not an option when the weight gets large.

The library syntax is `GEN `

.**zetamultinit**(long maxw, long prec)