Permutations are represented in gp as `t_VECSMALL`

s and can be input
directly as `Vecsmall([1,3,2,4])`

or obtained from the iterator
`forperm`

:

? forperm(3, p, print(p)) \\ iterate through S_{3}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 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 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 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} evaluated at a (`'x`

by default),
defined by

∑_{n = 0}^{ oo }B_{n}(x)(T^{n})/(n!) = (Te^{xT})/(e^{T}-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 B_{k}
at different values.

The library syntax is `GEN `

.
The variant **bernpol_eval**(long n, GEN a = NULL)`GEN `

returns
the k-the Bernoulli polynomial in variable v.**bernpol**(long k, long v)

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 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 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 = 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 `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 two
hours), the practical time complexity is closer to Õ(n^{log2 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,
lim_{t → 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} a^{n-k} b^{k}.
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 E_{n},
where E_{0} = 1, E_{1} = 0, E_{2} = -1,..., are integers such that
(1)/(cosh t) = ∑_{n ≥ 0} (E_{n})/(n!) t^{n}.
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 A_{n} in variable v defined by

∑_{n = 0}^{ oo } A_{n}(x) (T^{n})/(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 `

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

is a variable number.

Euler polynomial E_{n} in variable v defined by

∑_{n = 0}^{ oo } E_{n}(x)(T^{n})/(n!) = (2e^{xT})/(e^{T}+1).

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

The library syntax is `GEN `

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

is a variable number.

Euler number E_{n},
where E_{0} = 1, E_{1} = 0, E_{2} = -1,..., are integers such that
(1)/(cosh t) = ∑_{n ≥ 0} (E_{n})/(n!) t^{n}.
The argument n should be a nonnegative integer. Return E_{n}
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 E_{0},
E_{2},..., E_{2n}:

? eulervec(5) \\ E_{0}, E_{2}..., 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 `eulerfrac`

s
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
H_{n} = ∑_{i = 1}^{n} (1)/(i).
In general, this is
H_{n,r} = ∑_{i = 1}^{n} (1)/(i^{r}).
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 H_{n} as a
floating point number, use H_{n} = `psi`

(n+1) + `Euler`

.

The library syntax is `GEN `

.
Also available is **harmonic0**(ulong n, GEN r = NULL)`GEN `

for r = 1.**harmonic**(ulong n)

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 < 10^{15} (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([4]), 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) x^{k} = 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) x^{n} = (x^{k})/((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, B_{n}= 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 `

.
Also available are **stirling**(long n, long k, long flag)`GEN `

(**stirling1**(ulong n, ulong k)*flag* = 1) and `GEN `

(**stirling2**(ulong n, ulong k)*flag* = 2).