Pari/GP Reference Documentation Contents - Global index - GP keyboard shortcuts

## Transcendental functions

 Catalan   Euler   I   Pi   ^   abs   acos   acosh   agm   arg   asin   asinh   atan   atanh   bernfrac   bernpol   bernreal   bernvec   besselh1   besselh2   besseli   besselj   besseljh   besselk   besseln   cos   cosh   cotan   cotanh   dilog   eint1   erfc   eta   exp   expm1   gamma   gammah   gammamellininv   gammamellininvasymp   gammamellininvinit   hyperu   incgam   incgamc   lambertw   lngamma   log   polylog   psi   sin   sinc   sinh   sqr   sqrt   sqrtn   tan   tanh   teichmuller   theta   thetanullk   weber   zeta   zetamult   zetamultall   zetamultconvert 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, 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.10.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. If x is a p-adic number, its precision will increase if v_p(n) > 0. Powering a binary quadratic form (types `t_QFI` and `t_QFR`) returns a representative of the class, which is always reduced if the input was. (In particular, `x^1` returns x itself, whether it is reduced or not.) PARI is able to rewrite the multiplication x * x of two identical objects as x^2, or `sqr`(x). Here, identical means the operands are two different labels referencing the same chunk of memory; no equality test is performed. This is no longer true when more than two arguments are involved. * If the exponent n is not an integer, 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.) * As an exception, if the exponent is a rational number p/q and x an integer modulo a prime or a p-adic number, return a solution y of y^q = x^p if it exists. Currently, q must not have large prime factors. Beware that ``` ? Mod(7,19)^(1/2) %1 = Mod(11, 19) /* is any square root */ ? sqrt(Mod(7,19)) %2 = Mod(8, 19) /* is the smallest square root */ ? Mod(7,19)^(3/5) %3 = Mod(1, 19) ? %3^(5/3) %4 = Mod(1, 19) /* Mod(7,19) is just another cubic root */ ``` * If the exponent is a negative integer, an inverse must be computed. For non-invertible `t_INTMOD` x, this will fail and implicitly exhibit a non trivial 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 behaviour 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 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 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] ``` The library syntax is `GEN gpow(GEN x, GEN n, long prec)` for x^n. Catalan 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 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)`. I The complex number sqrt{-1}. The library syntax is `GEN gen_I()`. Pi 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)`. abs 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)`. acos 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)`. acosh 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)`. agm 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)`. arg Argument of the complex number x, such that -π < arg(x) ≤ π. The library syntax is `GEN garg(GEN x, long prec)`. asin 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)`. asinh 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)`. atan 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)`. atanh 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)`. bernfrac Bernoulli number B_x, where B_0 = 1, B_1 = -1/2, B_2 = 1/6,..., expressed as a rational number. The argument x should be of type integer. The library syntax is `GEN bernfrac(long x)`. bernpol 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 bernpol(long n, long v = -1)` where `v` is a variable number. bernreal Bernoulli number B_x, as `bernfrac`, but B_x is returned as a real number (with the current precision). The library syntax is `GEN bernreal(long x, long prec)`. bernvec This routine is obsolete, kept for backward compatibility only. The library syntax is `GEN bernvec(long x)`. besselh1 H^1-Bessel function of index nu and argument x. The library syntax is `GEN hbessel1(GEN nu, GEN x, long prec)`. besselh2 H^2-Bessel function of index nu and argument x. The library syntax is `GEN hbessel2(GEN nu, GEN x, long prec)`. besseli 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)`. besselj 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)`. besseljh 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.10.0, this function is not very accurate when x is small. The library syntax is `GEN jbesselh(GEN n, GEN x, long prec)`. besselk K-Bessel function of index nu and argument x. The library syntax is `GEN kbessel(GEN nu, GEN x, long prec)`. besseln N-Bessel function of index nu and argument x. The library syntax is `GEN nbessel(GEN nu, GEN x, long prec)`. cos Cosine of x. The library syntax is `GEN gcos(GEN x, long prec)`. cosh Hyperbolic cosine of x. The library syntax is `GEN gcosh(GEN x, long prec)`. cotan Cotangent of x. The library syntax is `GEN gcotan(GEN x, long prec)`. cotanh Hyperbolic cotangent of x. The library syntax is `GEN gcotanh(GEN x, long prec)`. dilog 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)`. eint1 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)`. erfc 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)`. eta 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). exp 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. expm1 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 = 28 ms. ? default(realprecision, 10040); x = 1e-100; ? c = expm1(x); \\ reference point ? abs(a-c)/c \\ relative error in expm1(x) %7 = 0.E-10017 ? 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 both faster and more accurate than `exp(x)-1`. The library syntax is `GEN gexpm1(GEN x, long prec)`. gamma 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 ggamma(GEN s, long prec)`. For a `t_PADIC` x, the function `GEN Qp_gamma(GEN x)` is also available. gammah Gamma function evaluated at the argument x+1/2. The library syntax is `GEN ggammah(GEN x, long prec)`. gammamellininv 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)`. gammamellininvasymp 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)`. gammamellininvinit 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)`. hyperu U-confluent hypergeometric function with parameters a and b. The parameters a and b can be complex but the present implementation requires x to be positive. The library syntax is `GEN hyperu(GEN a, GEN b, GEN x, long prec)`. incgam 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 incgam0(GEN s, GEN x, GEN g = NULL, long prec)`. Also available is `GEN incgam(GEN s, GEN x, long prec)`. incgamc 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)`. lambertw Lambert W function, solution of the implicit equation xe^x = y, for y > 0. The library syntax is `GEN glambertW(GEN y, long prec)`. lngamma 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)`. log 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. polylog 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 polylog0(long m, GEN x, long flag, long prec)`. Also available is `GEN gpolylog(long m, GEN x, long prec)` (flag = 0). psi The ψ-function of x, i.e. the logarithmic derivative Γ'(x)/Γ(x). The library syntax is `GEN gpsi(GEN x, long prec)`. sin Sine of x. The library syntax is `GEN gsin(GEN x, long prec)`. sinc 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)`. sinh Hyperbolic sine of x. The library syntax is `GEN gsinh(GEN x, long prec)`. sqr 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)`. sqrt 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. sqrtn 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. tan Tangent of x. The library syntax is `GEN gtan(GEN x, long prec)`. tanh Hyperbolic tangent of x. The library syntax is `GEN gtanh(GEN x, long prec)`. teichmuller 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^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)`. theta 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)`. thetanullk 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^iθ)/(dz^i)(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`). weber 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 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)`. zeta For s 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. For s 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). The library syntax is `GEN gzeta(GEN s, long prec)`. zetamult For s a vector of positive integers such that s[1] ≥ 2, returns the multiple zeta value (MZV) ζ(s_1,..., s_k) = ∑_{n_1 > ... > n_k > 0} n_1^{-s_1}...n_k^{-s_k}. ``` ? 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). In addition to the above format (`avec`), the function also accepts an internal binary format `evec` (each s_i as 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 zetamult(GEN s, long prec)`. zetamultall 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)`. zetamultconvert `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 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)`.