## Combinatorics

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 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`.

#### bernfrac(n)

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 `bervec` 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)`.

#### bernpol(n, {v = 'x})

Bernoulli polynomial Bn 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(n)

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)`.

#### bernvec(n)

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 `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 Õ(nlog2 6).

The library syntax is `GEN bernvec(long n)`.

#### binomial(n, {k})

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 Pochammer symbol.

```  ? 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 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)
%3 = [1, 4, 6, 4, 1]
```

The library syntax is `GEN binomial0(GEN n, GEN k = NULL)`.

#### eulerfrac(n)

Euler number En, where E0 = 1, E1 = 0, E2 = -1,..., are integers such that (1)/(cosh t) = ∑n ≥ 0 (En)/(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)`.

#### eulerianpol(n, {v = 'x})

Eulerian polynomial An in variable v.

```  ? 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.

#### eulerpol(n, {v = 'x})

Euler polynomial En in variable v.

```  ? 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.

#### eulerreal(n)

Euler number En, where E0 = 1, E1 = 0, E2 = -1,..., are integers such that (1)/(cosh t) = ∑n ≥ 0 (En)/(n!) t^n. 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)`.

#### eulervec(n)

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 `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)`.

#### fibonacci(x)

x-th Fibonacci number.

The library syntax is `GEN fibo(long x)`.

#### hammingweight(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)`.

#### harmonic(n, {r = 1})

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 = 1^n (1)/(i). In general, this is Hn,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 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.

#### numbpart(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 < 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)`.

#### numtoperm(n, k)

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)`.

#### partitions(k, {a = k}, {n = 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(), 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(), Vecsmall([1, 3]), Vecsmall([2, 2])]
? partitions(4,[0,3], 2) \\ at most 2 parts
%3 = [Vecsmall(), 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)`.

#### permcycles(x)

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(),Vecsmall(),Vecsmall()]
? permcycles(Vecsmall([2,3,1]))
%2 = [Vecsmall([1,2,3])]
? permcycles(Vecsmall([2,1,3]))
%3 = [Vecsmall([1,2]),Vecsmall()]
```

The library syntax is `GEN permcycles(GEN x)`.

#### permorder(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)`.

#### permsign(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)`.

#### permtonum(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(n, k, {flag = 1})

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; simplify(t));
}

/* Bell numbers, Bn = B[n+1] = sum(k = 0, n, S(n,k)), n = 0..N */
vecbell(N)=
{ my (B = vector(N+1));
B = B = 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).