Bernoulli number B_n, where B₀ = 1, B₁ = -1/2, B₂ = 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 B_n evaluated at a ('x by default), defined by

∑_{n = 0} ^oo B_n(x)(Tⁿ)/(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 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 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₀, B₂,..., B_2n:

  ? 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 Õ(n²); in the feasible range n ≤ 10⁵ (requires about two hours), the practical time complexity is closer to Õ(n^{log₂ 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)ⁿ = ∑_{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₀ = 1, E₁ = 0, E₂ = -1,..., are integers such that (1)/(cosh t) = ∑_{n ≥ 0} (E_n)/(n!) tⁿ. 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!) = (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 E_n in variable v defined by

∑_{n = 0} ^oo E_n(x)(Tⁿ)/(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 eulerpol(long n, long v = -1) where v is a variable number.

Euler number E_n, where E₀ = 1, E₁ = 0, E₂ = -1,..., are integers such that (1)/(cosh t) = ∑_{n ≥ 0} (E_n)/(n!) tⁿ. 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₀, E₂,..., E_2n:

  ? 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 H_n = ∑_{i = 1}ⁿ (1)/(i). In general, this is H_n,r = ∑_{i = 1}ⁿ (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 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 < 10¹⁵ (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) 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ⁿ = (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, 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).