Permutations are represented in gp as t_VECSMALLs and can be input
directly as Vecsmall([1,3,2,4]) or obtained from the iterator
forperm:
? forperm(3, p, print(p)) \\ iterate through S3 Vecsmall([1, 2, 3]) Vecsmall([1, 3, 2]) Vecsmall([2, 1, 3]) Vecsmall([2, 3, 1]) Vecsmall([3, 1, 2]) Vecsmall([3, 2, 1])
Permutations can be multiplied via *, raised to some power using
^, inverted using ^(-1), conjugated as
p * q * p^(-1). Their order and signature are available via
permorder and permsign.
Bernoulli number Bn,
where B0 = 1, B1 = -1/2, B2 = 1/6,..., expressed as a rational
number.
The argument n should be a nonnegative integer. The function bernvec
creates a cache of successive Bernoulli numbers which greatly speeds up
later calls to bernfrac:
? bernfrac(20000); time = 107 ms. ? bernvec(10000); \\ cache B0, B2, ..., B_20000 time = 35,957 ms. ? bernfrac(20000); \\ now instantaneous ?
The library syntax is GEN bernfrac(long n).
Bernoulli polynomial Bn evaluated at a ('x by default),
defined by
∑n = 0 oo Bn(x)(Tn)/(n!) = (TexT)/(eT-1).
? bernpol(1) %1 = x - 1/2 ? bernpol(3) %2 = x^3 - 3/2*x^2 + 1/2*x ? bernpol(3, 2) %3 = 3
Note that evaluation at a is only provided for convenience
and uniformity of interface: contrary to, e.g., polcyclo, computing
the evaluation is no faster than
B = bernpol(k); subst(B, 'x, a)
and the latter allows to reuse B to evaluate Bk at different values.
The library syntax is GEN bernpol_eval(long n, GEN a = NULL).
The variant GEN bernpol(long k, long v) returns
the k-the Bernoulli polynomial in variable v.
Bernoulli number
Bn, as bernfrac, but Bn is returned as a real number
(with the current precision). The argument n should be a nonnegative
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 B0, B2,..., B2n:
? bernvec(5) \\ B0, B2..., 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 = 18,245 ms. ? bernvec(5000); time = 1,338 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 bernfracs
after the bernvec call gives:
? forstep(n = 2, 10000, 2, bernfrac(n)) time = 1 ms.
The time and space complexity of this function are Õ(n2); in the feasible range n ≤ 105 (requires about two hours), the practical time complexity is closer to Õ(nlog2 6).
The library syntax is GEN bernvec(long n).
binomial coefficient binom{n}{k}. Here k must be an integer, but n can be any PARI object. For nonnegative k, binom{n}{k} = (n)k / k! is polynomial in n, where (n)k = n(n-1)...(n-k+1) is the Pochhammer symbol used by combinatorists (which is different from the one used by analysts).
? binomial(4,2)
%1 = 6
? n = 4; vector(n+1, k, binomial(n,k-1))
%2 = [1, 4, 6, 4, 1]
? binomial('x, 2)
%3 = 1/2*x^2 - 1/2*x
When n is a negative integer and k is negative,
we use Daniel Loeb's extension,
limt → 1 Γ(n+t) / Γ(k+t) / Γ(n-k+t).
(Sets with a negative number of elements, Adv. Math. 91 (1992),
no. 1, 64–74. See also https://arxiv.org/abs/1105.3689.)
This way the symmetry relation binom{n}{k} = binom{n}{n - k}
becomes valid for all integers n and k, and
the binomial theorem
holds for all complex numbers a, b, n with |b| < |a|:
(a + b)n = ∑k ≥ 0 binom{n}{k} an-k bk.
Beware that this extension is incompatible with another traditional
extension (binom{n}{k} := 0 if k < 0); to enforce the latter, use
BINOMIAL(n, k) = if (k >= 0, binomial(n, k));
The argument k may be omitted if n is a
nonnegative integer; in this case, return the vector with n+1
components whose k+1-th entry is binomial(n,k)
? binomial(4) %4 = [1, 4, 6, 4, 1]
The library syntax is GEN binomial0(GEN n, GEN k = NULL).
Euler number En, where E0 = 1, E1 = 0, E2 = -1,..., are integers such that (1)/(cosh t) = ∑n ≥ 0 (En)/(n!) tn. The argument n should be a nonnegative integer.
? vector(10,i,eulerfrac(i)) %1 = [0, -1, 0, 5, 0, -61, 0, 1385, 0, -50521] ? eulerfrac(20000); ? sizedigit(%)) %3 = 73416
The library syntax is GEN eulerfrac(long n).
Eulerian polynomial An in variable v defined by
∑n = 0 oo An(x) (Tn)/(n!) = (x-1)/(x-e(x-1)T).
? eulerianpol(2) %1 = x + 1 ? eulerianpol(5, 't) %2 = t^4 + 26*t^3 + 66*t^2 + 26*t + 1
The library syntax is GEN eulerianpol(long n, long v = -1) where v is a variable number.
Euler polynomial En in variable v defined by
∑n = 0 oo En(x)(Tn)/(n!) = (2exT)/(eT+1).
? eulerpol(1) %1 = x - 1/2 ? eulerpol(3) %2 = x^3 - 3/2*x^2 + 1/4
The library syntax is GEN eulerpol(long n, long v = -1) where v is a variable number.
Euler number En, where E0 = 1, E1 = 0, E2 = -1,..., are integers such that (1)/(cosh t) = ∑n ≥ 0 (En)/(n!) tn. The argument n should be a nonnegative integer. Return En as a real number (with the current precision).
? sizedigit(eulerfrac(20000)) %1 = 73416 ? eulerreal(20000); %2 = 9.2736664576330851823546169139003297830 E73414
The library syntax is GEN eulerreal(long n, long prec).
Returns a vector containing the nonzero Euler numbers E0, E2,..., E2n:
? eulervec(5) \\ E0, E2..., E_10 %1 = [1, -1, 5, -61, 1385, -50521] ? eulerfrac(10) %2 = -50521
This routine uses more memory but is faster than
repeated calls to eulerfrac:
? forstep(n = 2, 8000, 2, eulerfrac(n)) time = 27,3801ms. ? eulervec(4000); time = 8,430 ms.
The computed Euler numbers are stored in an incremental
cache which makes later calls to eulerfrac and eulerreal
instantaneous in the cache range: re-running the same previous eulerfracs
after the eulervec call gives:
? forstep(n = 2, 10000, 2, eulerfrac(n)) time = 0 ms.
The library syntax is GEN eulervec(long n).
x-th Fibonacci number.
The library syntax is GEN fibo(long x).
If x is a t_INT, return the binary Hamming weight of |x|. Otherwise
x must be of type t_POL, t_VEC, t_COL, t_VECSMALL, or
t_MAT and the function returns the number of nonzero coefficients of
x.
? hammingweight(15) %1 = 4 ? hammingweight(x^100 + 2*x + 1) %2 = 3 ? hammingweight([Mod(1,2), 2, Mod(0,3)]) %3 = 2 ? hammingweight(matid(100)) %4 = 100
The library syntax is long hammingweight(GEN x).
Generalized harmonic number of index n ≥ 0 in power r, as a rational number. If r = 1 (or omitted), this is the harmonic number Hn = ∑i = 1n (1)/(i). In general, this is Hn,r = ∑i = 1n (1)/(ir). The function runs in time Õ(r n), essentially linear in the size of the output.
? harmonic(0) %1 = 0 ? harmonic(1) %2 = 1 ? harmonic(10) %3 = 7381/2520 ? harmonic(10, 2) %4 = 1968329/1270080 ? harmonic(10, -2) %5 = 385
Note that the numerator and denominator are of order
exp((r+o(1))n) and this will overflow for large n. To obtain Hn as a
floating point number, use Hn = psi(n+1) + Euler.
The library syntax is GEN harmonic0(ulong n, GEN r = NULL).
Also available is GEN harmonic(ulong n) for r = 1.
Gives the number of unrestricted partitions of
n, usually called p(n) in the literature; in other words the number of
nonnegative integer solutions to a+2b+3c+.. .= n. n must be of type
integer and n < 1015 (with trivial values p(n) = 0 for n < 0 and
p(0) = 1). The algorithm uses the Hardy-Ramanujan-Rademacher formula.
To explicitly enumerate them, see partitions.
The library syntax is GEN numbpart(GEN n).
Generates the k-th permutation (as a row vector of length n) of the
numbers 1 to n. The number k is taken modulo n!, i.e. inverse
function of permtonum. The numbering used is the standard lexicographic
ordering, starting at 0.
The library syntax is GEN numtoperm(long n, GEN k).
Returns the vector of partitions of the integer k as a sum of positive
integers (parts); for k < 0, it returns the empty set [], and for k
= 0 the trivial partition (no parts). A partition is given by a
t_VECSMALL, where parts are sorted in nondecreasing order:
? partitions(3) %1 = [Vecsmall([3]), Vecsmall([1, 2]), Vecsmall([1, 1, 1])]
correspond to 3, 1+2 and 1+1+1. The number
of (unrestricted) partitions of k is given
by numbpart:
? #partitions(50) %1 = 204226 ? numbpart(50) %2 = 204226
Optional parameters n and a are as follows:
* n = nmax (resp. n = [nmin,nmax]) restricts partitions to length less than nmax (resp. length between nmin and nmax), where the length is the number of nonzero entries.
* a = amax (resp. a = [amin,amax]) restricts the parts to integers less than amax (resp. between amin and amax).
? partitions(4, 2) \\ parts bounded by 2 %1 = [Vecsmall([2, 2]), Vecsmall([1, 1, 2]), Vecsmall([1, 1, 1, 1])] ? partitions(4,, 2) \\ at most 2 parts %2 = [Vecsmall([4]), Vecsmall([1, 3]), Vecsmall([2, 2])] ? partitions(4,[0,3], 2) \\ at most 2 parts %3 = [Vecsmall([1,3]), Vecsmall([2,2])]
By default, parts are positive and we remove zero entries unless amin ≤ 0, in which case nmin is ignored and we fix #X = nmax:
? partitions(4, [0,3]) \\ parts between 0 and 3
%1 = [Vecsmall([0, 0, 1, 3]), Vecsmall([0, 0, 2, 2]),\
Vecsmall([0, 1, 1, 2]), Vecsmall([1, 1, 1, 1])]
? partitions(1, [0,3], [2,4]) \\ no partition with 2 to 4 nonzero parts
%2 = []
The library syntax is GEN partitions(long k, GEN a = NULL, GEN n = NULL).
Given a permutation x on n elements, return the orbits of {1,...,n} under the action of x as cycles.
? permcycles(Vecsmall([1,2,3])) %1 = [Vecsmall([1]),Vecsmall([2]),Vecsmall([3])] ? permcycles(Vecsmall([2,3,1])) %2 = [Vecsmall([1,2,3])] ? permcycles(Vecsmall([2,1,3])) %3 = [Vecsmall([1,2]),Vecsmall([3])]
The library syntax is GEN permcycles(GEN x).
Given a permutation x on n elements, return its order.
? p = Vecsmall([3,1,4,2,5]); ? p^2 %2 = Vecsmall([4,3,2,1,5]) ? p^4 %3 = Vecsmall([1,2,3,4,5]) ? permorder(p) %4 = 4
The library syntax is GEN permorder(GEN x).
Given a permutation x on n elements, return its signature.
? p = Vecsmall([3,1,4,2,5]); ? permsign(p) %2 = -1 ? permsign(p^2) %3 = 1
The library syntax is long permsign(GEN x).
Given a permutation x on n elements, gives the number k such that
x = numtoperm(n,k), i.e. inverse function of numtoperm.
The numbering used is the standard lexicographic ordering, starting at 0.
The library syntax is GEN permtonum(GEN x).
Stirling number of the first kind s(n,k) (flag = 1, default) or of the second kind S(n,k) (flag = 2), where n, k are nonnegative integers. The former is (-1)n-k times the number of permutations of n symbols with exactly k cycles; the latter is the number of ways of partitioning a set of n elements into k nonempty subsets. Note that if all s(n,k) are needed, it is much faster to compute ∑k s(n,k) xk = x(x-1)...(x-n+1). Similarly, if a large number of S(n,k) are needed for the same k, one should use ∑n S(n,k) xn = (xk)/((1-x)...(1-kx)). (Should be implemented using a divide and conquer product.) Here are simple variants for n fixed:
/* list of s(n,k), k = 1..n */
vecstirling(n) = Vec( factorback(vector(n-1,i,1-i*'x)) )
/* list of S(n,k), k = 1..n */
vecstirling2(n) =
{ my(Q = x^(n-1), t);
vector(n, i, t = divrem(Q, x-i); Q=t[1]; simplify(t[2]));
}
/* Bell numbers, Bn = B[n+1] = sum(k = 0, n, S(n,k)), n = 0..N */
vecbell(N)=
{ my (B = vector(N+1));
B[1] = B[2] = 1;
for (n = 2, N,
my (C = binomial(n-1));
B[n+1] = sum(k = 1, n, C[k]*B[k]);
); B;
}
The library syntax is GEN stirling(long n, long k, long flag).
Also available are GEN stirling1(ulong n, ulong k)
(flag = 1) and GEN stirling2(ulong n, ulong k) (flag = 2).