Catalan's constant G = ∑_{n >= 0}((-1)ⁿ)/((2n+1)²) = 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 genI().

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²}). 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²}). 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²}). 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).

H¹-Bessel function of index nu and argument x.

The library syntax is GEN hbessel1(GEN nu, GEN x, long prec).

H²-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.17.0, this function is not very accurate when x is small.

The library syntax is GEN jbesselh(GEN n, GEN x, long prec).

k-th zero of the J-Bessel function of index nu, close to π(ν/2 + k - 1/4), usually noted j_ν,k.

  ? besseljzero(0) \\  {first zero of J0}
  %1 = 2.4048255576957727686216318793264546431
  ? besselj(0, %)
  %2 = 7.1951595399463653939930598011247182898 E-41
  ? besseljzero(0, 2) \\  {second zero}
  %3 = 5.5200781102863106495966041128130274252
  ? besseljzero(I) \\  {also works for complex orders, here Ji}
  %4 = 2.5377... + 1.4753...*I

The function uses a Newton iteration due to Temme. If ν is real and nonnegative, the result is guaranteed and the function really returns the k-th positive zero of J_ν. For general ν, the result is not well defined, given by the Newton iteration with π(ν/2 + k - 1/4) as a starting value. (N.B. Using this method for large real ν would give completely different results than the j_ν,k unless k is large enough.)

The library syntax is GEN besseljzero(GEN nu, long k, long bitprec).

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 ybessel(GEN nu, GEN x, long prec).

Y-Bessel function of index nu and argument x.

The library syntax is GEN ybessel(GEN nu, GEN x, long prec).

k-th zero of the Y-Bessel function of index nu, close to π(ν/2 + k - 3/4), usually noted y_ν,k.

  ? besselyzero(0) \\  {first zero of Y0}
  %1 = 0.89357696627916752158488710205833824123
  ? bessely(0, %)
  %2 = 1.8708573650996561952 E-39
  ? besselyzero(0, 2) \\  {second zero}
  %3 = 3.9576784193148578683756771869174012814
  ? besselyzero(I) \\  {also works for complex orders, here Yi}
  %4 = 1.03930... + 1.3266...*I

The function uses a Newton iteration due to Temme. If ν is real and nonnegative, the result is guaranteed and the function really returns the k-th positive zero of Y_ν. For general ν, the result is not well defined, given by Newton iteration with π(ν/2 + k - 3/4) as a starting value. (N.B. Using this method for large real ν would give completely different results than the y_ν,k unless k is large enough.)

The library syntax is GEN besselyzero(GEN nu, long k, long bitprec).

Cosine of x. Note that, for real x, cosine and sine can be obtained simultaneously as

  cs(x) = my(z = exp(I*x)); [real(z), imag(z)];

and for general complex x as

  cs2(x) = my(z = exp(I*x), u = 1/z); [(z+u)/2, (z-u)/2];

Note that the latter function suffers from catastrophic cancellation when z² ~ ±1.

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 Li₂(x) = ∑_{n ≥ 1}xⁿ/n².

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 veceint1(GEN x, GEN n = NULL, long prec). Also available is GEN eint1(GEN x, long prec).

Complete elliptic integral of the second kind E(k) = ∫₀^π/2(1-k²sin(t)²)^1/2dt for the complex parameter k using the agm.

In particular, the perimeter of an ellipse of semi-major and semi-minor axes a and b is given by

    e = sqrt(1 - (b/a)^2); \\ eccentricity
    4 * a * ellE(e)  \\ perimeter

The library syntax is GEN ellE(GEN k, long prec).

Complete elliptic integral of the first kind K(k) = ∫₀^π/2(1-k²sin(t)²)^-1/2dt for the complex parameter k using the agm.

The library syntax is GEN ellK(GEN k, long prec).

Complementary error function, analytic continuation of (2/sqrtπ)∫_x ^oo e^{-t^{2}}dt = {sign(x)}incgam(1/2,x²)/sqrtπ for real x ! = 0. The latter expression extends the function definition from real x to complex x with positive real part (or zero real part and positive imaginary part). This is extended to the whole complex plane by the functional equation erfc(-x) = 2 - erfc(x).

  ? erfc(0)
  %1 = 1.0000000000000000000000000000000000000
  ? erfc(1)
  %2 = 0.15729920705028513065877936491739074071
  ? erfc(1+I)
  %3 = -0.31615128169794764488027108024367036903
       - 0.19045346923783468628410886196916244244*I

The library syntax is GEN gerfc(GEN x, long prec).

Variants of Dedekind's η function. If flag = 0, return ∏_{n = 1} ^oo (1-qⁿ), 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 nonzero, x is converted to a complex number and we return the true η function, q^1/24∏_{n = 1} ^oo (1-qⁿ), 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 gexp(GEN x, long prec). For a t_PADIC x, the function GEN Qp_exp(GEN x) is also available.

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, which is the analytic continuation of Γ(s) = ∫₀ ^oo t^s-1exp(-t)dt, Re(s) > 0. Error if s is a nonpositive integer, where Γ has a (simple) pole.

  ? gamma(5)  \\   Γ(n) = (n-1)! for a positive integer n
  %1 = 24.000000000000000000000000000000000000
  ? gamma(0)
   ***   at top-level: gamma(0)
   ***                 ^ —  — --
   *** gamma: domain error in gamma: argument = nonpositive integer
  
  ? gamma(x + O(x^3))
  %2 = x^-1 - 0.57721566490153286060651209008240243104 + O(x)

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 ggamma(GEN s, long prec). For a t_PADIC x, the function GEN Qp_gamma(GEN x) is also available.

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. If the optional parameter m is present, return the m-th derivative G^(m)(t).

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

The 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₁)...Γ_ℝ(s+a_d) , where A is the vector [a₁,...,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. We also allow A to be the output of gammamellininvinit.

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₁)...Γ_ℝ(s+a_d) where A is the vector [a₁,...,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. By extension, A is allowed to be an Ldata or an Linit, understood as the inverse Mellin transform of the L-function gamma factor.

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²):

  ? 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²), and its second derivative is 4π z exp(-π z²)(2π z² - 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 argument z, which may be complex, p-adic or a power series.

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!) , where (a)_n = a(a+1)...(a+n-1) is the rising Pochhammer 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 ₀F₁(;b;z).

The non-archimedean cases (z a p-adic or power series) are handled trivially. We now discuss the case of a complex 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, by analytic continuation.

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 .

F10: hypergeom(a,,z) This is (1-z)^-a.

F11: hypergeom(a,b,z) is the Kummer confluent hypergeometric function, computed by summing the series. R = oo

F20: hypergeom([a,b],,z). R = 0, computed as (1)/(Γ(a))∫₀ ^oo t^a-1(1-zt)^-be^-tdt .

F21: hypergeom([a,b],c,z) (or [c]). R = 1, extended by (Γ(c))/(Γ(b)Γ(c-b)) ∫₀¹ 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.

F31: hypergeom([a,b,c],d,z) (or [d]). R = 0, computed as (1)/(Γ(a))∫₀ ^oo t^a-1e^-t ₂F₁(b,c;d;tz)dt .

F32: hypergeom([a,b,c],[d,e],z). R = 1, extended by (Γ(e))/(Γ(c)Γ(e-c)) ∫₀¹t^c-1(1-t)^e-c-1₂F₁(a,b;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 ₂F₀(a,b,z) = (-z)^-aU(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^-tt^s-1dt, 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 incgam0(GEN s, GEN x, GEN g = NULL, long prec). Also available is 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 ∫₀^× e^-tt^s-1dt.

The library syntax is GEN incgamc(GEN s, GEN x, long prec).

Lambert W function, solution of the implicit equation xe^× = y.

* For real inputs y: If branch = 0, principal branch W₀ defined for y ≥ -exp(-1). If branch = -1, branch W_-1 defined for -exp(-1) ≤ y < 0.

* For p-adic inputs, p odd: give a solution of xexp(x) = y if y has positive valuation, of log(x)+x = log(y) otherwise.

* For 2-adic inputs: give a solution of xexp(x) = y if y has valuation > 1, of log(x)+x = log(y) otherwise.

Caveat. Complex values of y are also supported but experimental. The other branches W_k for k not equal to 0 or -1 (set branch to k) are also experimental.

For k ≥ 1, W_-1-k(x) = W_k(x), and Im(W_k(x)) is close to (π/2)(4k-sign(x)).

The library syntax is GEN glambertW(GEN y, long branch, long prec).

Lerch transcendent Φ(z,s,a) = ∑_{n ≥ 0}zⁿ(n+a)^-s and analytically continued, for reasonable values of the arguments.

The library syntax is GEN lerchphi(GEN z, GEN s, GEN a, long prec).

Lerch zeta function L(s,a,λ) = ∑_{n ≥ 0}e^{2π iλ n}(n+a)^-s and analytically continued, for reasonable values of the arguments.

The library syntax is GEN lerchzeta(GEN s, GEN a, GEN lam, long prec).

Principal branch of the logarithm of the gamma function of x. This function is analytic on the complex plane with nonpositive 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

For x a t_PADIC, return the p-adic logΓ_p function, which is the p-adic logarithm of Morita's gamma function for x ∈ ℤ_p, and Diamond's function if |x| > 1.

  ? lngamma(5+O(5^7))
  %2 = 4*5^2 + 4*5^3 + 5^4 + 2*5^5 + O(5^6)
  ? log(gamma(5+O(5^7)))
  %3 = 4*5^2 + 4*5^3 + 5^4 + 2*5^5 + O(5^6)
  ? lngamma(1/5+O(5^4))
  %4 = 4*5^-1 + 4 + 2*5 + 5^2 + 5^3 + O(5^4)
  ? gamma(1/5+O(5^4))
   ***   at top-level: gamma(1/5+O(5^4))
   ***                 ^ —  —  —  —  — --
   *** gamma: domain error in gamma: vp(x) < 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 glog(GEN x, long prec). For a t_PADIC x, the function GEN Qp_log(GEN x) is also available.

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^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|² ≤ 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-1B_m)/(m!)(log|x|)^m).

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

The library syntax is GEN polylog0(long m, GEN x, long flag, long prec). Also available is GEN gpolylog(long m, GEN x, long prec) (flag = 0).

For s a vector of positive integers and z a vector of complex numbers of the same length, returns the multiple polylogarithm value (MPV) ζ(s₁,..., s_r; z₁,...,z_r) = ∑_{n_{1} > ... > n_r > 0} ∏_{1 ≤ i ≤ r}z_i^n_i/n_i^s_i. If z is omitted, assume z = [1,...,1], i.e., Multiple Zeta Value. More generally, return Yamamoto's interpolation between ordinary multiple polylogarithms (t = 0) and star polylogarithms (t = 1, using the condition n₁ ≥ ... ≥ n_r > 0), evaluated at t.

We must have |z₁...z_i| ≤ 1 for all i, and if s₁ = 1 we must have z₁ ! = 1.

  ? 8*polylogmult([2,1],[-1,1]) - zeta(3)
  %1 = 0.E-38

Warning. The algorithm used converges when the z_i are ± 1. It may not converge as some z_i ! = 1 becomes too close to 1, even at roots of 1 of moderate order:

  ? polylogmult([2,1], (99+20*I)/101 * [1,1])
   *** polylogmult: sorry, polylogmult in this range is not yet implemented.
  ? polylogmult([2,1], exp(I*Pi/20)* [1,1])
   *** polylogmult: sorry, polylogmult in this range is not yet implemented.

More precisely, if y_i := 1 / (z₁...z_i) and v := min_{i < j; y_{i} ! = 1} |(1 - y_i) y_j| > 1/4 then the algorithm computes the value up to a 2^-b absolute error in O(k²N) operations on floating point numbers of O(N) bits, where k = ∑_i s_i is the weight and N = b / log₂ (4v).

The library syntax is GEN polylogmult_interpolate(GEN s, GEN z = NULL, GEN t = NULL, long prec). Also available is GEN polylogmult(GEN s, GEN z, long prec) (t is NULL).

The ψ-function of x, i.e. the logarithmic derivative Γ'(x)/Γ(x). If der is set, return the der-th derivative. For s a t_PADIC, evaluates the der-th derivative of the Morita ψ function at s.

The library syntax is GEN gpsi_der(GEN x, long der, long prec). For a t_PADIC x, the function GEN Qp_psi(GEN x, long der) is also available. For der = 0, GEN gpsi(GEN x, long prec) is also available.

Return the column vector v of all complex N-th roots of 1, where N is a positive integer. In other words, v[k] = exp(2I(k-1)π/N) for k = 1,..., N. Rational components (e.g., the roots ±1 and ± I) are given exactly, not as floating point numbers:

  ? rootsof1(4)
  %1 = [1, I, -1, -I]~
  ? rootsof1(3)
  %2 = [1, -1/2 + 0.866025...*I, -1/2 - 0.866025...*I]~

The library syntax is GEN grootsof1(long N, long prec).

Sine of x. Note that, for real x, cosine and sine can be obtained simultaneously as

  cs(x) = my(z = exp(I*x)); [real(z), imag(z)];

and for general complex x as

  cs2(x) = my(z = exp(I*x), u = 1/z); [(z+u)/2, (z-u)/2];

Note that the latter function suffers from catastrophic cancellation when z² ~ ±1.

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² = sinc(x/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 gsqrt(GEN x, long prec). For a t_PADIC x, the function GEN Qp_sqrt(GEN x) is also available.

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 gsqrtn(GEN x, GEN n, GEN *z = NULL, long prec). If x is a t_PADIC, the function GEN Qp_sqrtn(GEN x, GEN n, GEN *z) is also available.

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^v_p(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ⁿ) 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₂(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 precomputation 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 θ₁(z, q) = 2q^1/4 ∑_{n ≥ 0} (-1)ⁿ qⁿ⁽ⁿ⁺¹⁾ 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 vecthetanullk(GEN q, long k, long prec) returns the vector of all (dⁱθ)/(dzⁱ)(q,0) for all odd i = 1, 3,..., 2k-1. GEN vecthetanullk_tau(GEN tau, long k, long prec) returns 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²⁴-16)³/f²⁴, where j is the elliptic j-invariant (see the function ellj). If flag = 1, returns f₁(x) = η(x/2)/η(x) {such that} j = (f₁²⁴+16)³/f₁²⁴. Finally, if flag = 2, returns f₂(x) = sqrt{2}η(2x)/η(x) {such that} j = (f₂²⁴+16)³/f₂²⁴. Note the identities f⁸ = f₁⁸+f₂⁸ and ff₁f₂ = sqrt2.

The library syntax is GEN weber0(GEN x, long flag, long prec). Also available are GEN weberf(GEN x, long prec), GEN weberf1(GEN x, long prec) and 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 derivative can be used to compute Barnes' multiple gamma functions. For instance:

  ? mygamma(z)=exp(zetahurwitz(0,z,1)-zeta'(0));
  /* Alternate way to compute the gamma function */
  ? BarnesG(z)=exp(-zetahurwitz(-1,z,1)+(z-1)*lngamma(z)+zeta'(-1));
  /* Barnes G function, satisfying G(z+1)=gamma(z)*G(z): */
  ? BarnesG(6)/BarnesG(5)
  % = 24.000000000000000000000000000000000002

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₁,..., s_k) = ∑_{n_{1} > ... > n_k > 0} n₁^-s₁...n_k^-s_k of length k and weight ∑_i s_i. More generally, return Yamamoto's t-MZV interpolation evaluated at t: for t = 0, this is the ordinary MZV; for t = 1, we obtain the MZSV star value, with ≥ instead of strict inequalities; and of course, for t = 'x we obtain Yamamoto's one-variable polynomial.

  ? zetamult([2,1]) - zeta(3) \\ Euler's identity
  %1 = 0.E-38
  ? zetamult([2,1], 1)   \\ star value
  %2 = 2.4041138063191885707994763230228999815
  ? zetamult([2,1], 'x)
  %3 = 1.20205[...]*x + 1.20205[...]

If the bit precision is B, this function runs in time Õ(k(B+k)²) if t = 0, and Õ(kB³) otherwise.

In addition to the above format (avec), the function also accepts a binary word format evec (each s_i is replaced by s_i bits, all of them 0 but the last one) giving the MZV representation as an iterated integral, 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 zetamult_interpolate(GEN s, GEN t = NULL, long prec). Also available is GEN zetamult(GEN s, long prec) for t = 0.

List of all multiple zeta values (MZVs) for weight s₁ +...+ s_r up to k. Binary digits of flag mean : 0 = star values if set; 1 = values up to to duality if set (see zetamultdual, ignored if star values); 2 = values of weight k if set (else all values up to weight k); 3 = return the 2-component vector [Z, M], where M is the vector of the corresponding indices m, i.e., such that zetamult(M[i]) = Z[i]. Note that it is necessary to use zetamultconvert to have the corresponding avec (s₁,..., s_r).

With the default value flag = 0, the function returns a vector with 2^k-1-1 components whose i-th entry is the MZV of index i (see zetamult). If the bit precision is B, this function runs in time O(2^k k B²) for an output of size O(2^k B).

  ? Z = zetamultall(5); #Z \\ 2^4 - 1 MZVs of weight <= 5
  %1 = 15
  ? Z[10]
  %2 = 0.22881039760335375976874614894168879193
  ? zetamultconvert(10)
  %3 = Vecsmall([3, 2]) \\  {index 10 corresponds to ζ(3,2)}
  ? zetamult(%)  \\ double check
  %4 = 0.22881039760335375976874614894168879193
  ? zetamult(10) \\ we can use the index directly
  %5 = 0.22881039760335375976874614894168879193

If we use flag bits 1 and 2, we avoid unnecessary computations and copying, saving a potential factor 4: half the values are in lower weight and computing up to duality save another rough factor 2. Unfortunately, the indexing now no longer corresponds to the new shorter vector of MZVs:

  ? Z = zetamultall(5, 2); #Z \\ up to duality
  %6 = 9
  ? Z = zetamultall(5, 2); #Z \\ only weight 5
  %7 = 8
  ? Z = zetamultall(5, 2 + 4); #Z \\ both
  %8 = 4

So how to recover the value attached to index 10 ? Flag bit 3 returns the actual indices used:

  ? [Z, M] = zetamultall(5, 2 + 8); M \\ other indices were not included
  %9 = Vecsmall([1, 2, 4, 5, 6, 8, 9, 10, 12])
  ? Z[8] \\ index m = 10 is now in M[8]
  %10 = 0.22881039760335375976874614894168879193
  ? [Z, M] = zetamultall(5, 2 + 4 + 8); M
  %11 = Vecsmall([8, 9, 10, 12])
  ? Z[3] \\ index m = 10 is now in M[3]
  %12 = 0.22881039760335375976874614894168879193

The following construction automates the above programmatically, looking up the MZVs of index 10 ( = ζ(3,2)) in all cases, without inspecting the various index sets M visually:

  ? Z[vecsearch(M, 10)] \\ works in all the above settings
  %13 = 0.22881039760335375976874614894168879193

The library syntax is GEN zetamultall(long k, long flag, long prec).

a being either an evec, avec, or index m, converts into evec (flag = 0), avec (flag = 1), or index m (flag = 2).

  ? 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 flag).

s being either an evec, avec, or index m, return the dual sequence in avec format. The dual of a sequence of length r and weight k has length k-r and weight k. Duality is an involution and zeta values attached to dual sequences are the same:

  ? zetamultdual([4])
  %1 = Vecsmall([2, 1, 1])
  ? zetamultdual(%)
  %2 = Vecsmall([4])
  ? zetamult(%1) - zetamult(%2)
  %3 = 0.E-38

In evec form, duality simply reverses the order of bits and swaps 0 and 1:

  ? zetamultconvert([4], 0)
  %4 = Vecsmall([0, 0, 0, 1])
  ? zetamultconvert([2,1,1], 0)
  %5 = Vecsmall([0, 1, 1, 1])

The library syntax is GEN zetamultdual(GEN s).