Normal forms for products of gamma functions

This module rewrites an expression that is a finite product and quotient of gamma functions, factorials and binomials with rational-linear arguments a*n + b into a normal form built from factorials:

C * d^n * prod (A_i*n + B_i)!^e_i * prod p_k(n)^eps_k * prod gamma(a_j*n + b_j)^delta_j

where C is a constant (a rational times rational powers of primes and of pi, times gamma values at rational points in irreducible position), d is a nonzero rational, the factorial arguments are nonnegative integers for all n >= 0, and the p_k(n) are products of consecutive linear factors. The remaining gamma factors, if any, are in irreducible position: their shifts have pairwise distinct fractional parts within each slope, are normalized into the interval \([0, 1)\), and no complete “circle” is present.

For example:

sage: from sage.symbolic.gamma_products import simplify_gamma
sage: n = var('n')
sage: simplify_gamma(gamma(n + 2/3)*gamma(n + 1/3)/gamma(2/3)/gamma(1/3), n)
(1/27)^n*factorial(3*n)/factorial(n)

>>> from sage.all import *
>>> from sage.symbolic.gamma_products import simplify_gamma
>>> n = var('n')
>>> simplify_gamma(gamma(n + Integer(2)/Integer(3))*gamma(n + Integer(1)/Integer(3))/gamma(Integer(2)/Integer(3))/gamma(Integer(1)/Integer(3)), n)
(1/27)^n*factorial(3*n)/factorial(n)

from sage.symbolic.gamma_products import simplify_gamma
n = var('n')
simplify_gamma(gamma(n + 2/3)*gamma(n + 1/3)/gamma(2/3)/gamma(1/3), n)

A term can be written without gamma functions if and only if its fractional gamma shifts decompose into complete circles; the algorithm (see simplify_gamma()) searches for such a decomposition, and returns the remaining gamma factors in irreducible position otherwise. The search is subject to a budget; if it is exhausted, the returned normal form is irreducible only in the sense that no further circle applies, and it is not established that no other sequence of circle moves could eliminate the remaining gamma factors.

The same machinery evaluates finite products with a symbolic number of factors exactly, see factorial_product():

sage: k = var('k')
sage: from sage.symbolic.gamma_products import factorial_product
sage: f = (54*k^3 + 81*k^2 + 39*k + 6)/(2*k^3 + 9*k^2 + 13*k + 6)
sage: factorial_product(f, k, 0, n - 1)
factorial(3*n)/((2*n + 1)*(n + 1)*factorial(n)^3)

>>> from sage.all import *
>>> k = var('k')
>>> from sage.symbolic.gamma_products import factorial_product
>>> f = (Integer(54)*k**Integer(3) + Integer(81)*k**Integer(2) + Integer(39)*k + Integer(6))/(Integer(2)*k**Integer(3) + Integer(9)*k**Integer(2) + Integer(13)*k + Integer(6))
>>> factorial_product(f, k, Integer(0), n - Integer(1))
factorial(3*n)/((2*n + 1)*(n + 1)*factorial(n)^3)

k = var('k')
from sage.symbolic.gamma_products import factorial_product
f = (54*k^3 + 81*k^2 + 39*k + 6)/(2*k^3 + 9*k^2 + 13*k + 6)
factorial_product(f, k, 0, n - 1)

Products of more general bodies, whose rational part is possibly hidden, for example in ratios of factorials, are evaluated by exact_product(), which is used by symbolic_product():

sage: from sage.symbolic.gamma_products import exact_product
sage: exact_product(factorial(k)/factorial(k + 2)*2^k, k, 0, n)
2^(1/2*n^2 + 1/2*n)/((n + 2)*factorial(n + 1)^2)

>>> from sage.all import *
>>> from sage.symbolic.gamma_products import exact_product
>>> exact_product(factorial(k)/factorial(k + Integer(2))*Integer(2)**k, k, Integer(0), n)
2^(1/2*n^2 + 1/2*n)/((n + 2)*factorial(n + 1)^2)

from sage.symbolic.gamma_products import exact_product
exact_product(factorial(k)/factorial(k + 2)*2^k, k, 0, n)

ALGORITHM:

The expression E is hypergeometric, so the ratio E(n+1)/E(n) reduces, by applying the recurrence of the gamma function, to a rational function. Factoring the ratio into linear factors (A*n + B)^e, each factor contributes A^n * gamma(n + B/A)^e / gamma(B/A)^e, so everything becomes a finite multiset of gammas with slope 1 and rational shift, plus constants.

Fractional shifts are then eliminated with the Gauss multiplication formula: a complete circle {gamma(a*n + b + j/m) : j = 0..m-1} collapses to gamma(m*a*n + m*b), where members with integer argument are factorials and may be introduced freely. Every circle modulus divides the common denominator of the shifts present, and each move strictly decreases the total number of fractional gamma factors, counted with multiplicity, so the depth-first search with memoization terminates. The search is bounded; when no factorial form is found within the bound, or none exists, a normal form is committed instead: a state to which no circle applies.

The reduction of the term ratio to a rational function and its factorization into linear factors is the starting point of Gosper’s algorithm [Gos1978]; see also [PWZ1997], Chapter 5, and [GG1999], Section 23.4. The products of consecutive linear factors arising in the normal form correspond to greatest factorial factorizations [Paule1995], see also [GG1999], Section 23.3. That the three classical formulas generate all identities between products of gamma functions is folklore.

AUTHORS:

  • Martin Rubey (2026-10)

sage.symbolic.gamma_products.exact_product(f, var, a, b)[source]

Evaluate a finite product of a symbolic body exactly, as far as possible.

INPUT:

  • f – a symbolic expression, the body of the product

  • var – the product variable

  • a – the lower bound of the product

  • b – the upper bound of the product, such that the number of factors b - a + 1 is a nonnegative integer, or an integer-linear expression in a single variable

The body is first normalized to expose hidden rationality: factorials are rewritten as gamma functions, ratios of gamma functions whose arguments differ by integers are converted to polynomials, and binomials as well as falling and rising factorials with constant second argument are expanded. For example, the body factorial(k)/factorial(k+2) is rational in disguise and normalizes to 1/((k+1)*(k+2)).

The normalized body is then split multiplicatively into

  • a rational function of the product variable, which is evaluated in the factorial normal form described in simplify_gamma(),

  • powers c^p with base c independent of the product variable and polynomial exponent p, which are evaluated as c^symbolic_sum(p, var, a, b); for instance, 2^k over k = 0..n evaluates to 2^(n*(n+1)/2),

  • linear factors A*var + B with rational slope A and B free of the product variable, which are evaluated in terms of gamma functions, for instance x + k over k = 0..n evaluates to gamma(x+n+1)/gamma(x); over a concrete range the finite product is evaluated directly, since the gamma quotient has removable poles at the zeros of the factors, so that its size grows linearly with the number of factors,

  • factors independent of the product variable, which are evaluated as c^(b - a + 1),

  • a residual part, which is returned as an unevaluated product in the original variable and bounds.

The number of factors b - a + 1 must be a nonnegative integer, or an integer-linear expression q*n + r in a single variable n with q a positive integer and r a nonnegative integer; a ValueError is raised otherwise, and also if nothing can be evaluated. In both cases the caller may want to fall back to another method.

EXAMPLES:

sage: from sage.symbolic.gamma_products import exact_product
sage: k, n = var('k n')
sage: exact_product(2^k, k, 0, n)
2^(1/2*n^2 + 1/2*n)
sage: exact_product(factorial(k)/factorial(k + 2), k, 0, n)
1/((n + 2)*factorial(n + 1)^2)
sage: exact_product((3*k + 1)*(3*k + 2)*factorial(k), k, 0, n)
(1/3)^n*(3*n + 2)*(3*n + 1)*factorial(3*n)*product(factorial(k), k, 0, n)/factorial(n)

Powers with symbolic base and polynomial exponent, as well as
factors independent of the product variable, are evaluated; a zero
base is not aggregated, but evaluated at every position of a
concrete range, since it may have factors that are not defined::

    sage: x = var('x')
    sage: exact_product(x^k*(k + 1), k, 0, n)
    x^(1/2*n^2 + 1/2*n)*factorial(n + 1)
    sage: exact_product(sqrt(2)^k, k, 0, n)
    2^(1/4*n^2 + 1/4*n)
    sage: exact_product(0^k, k, 0, 3)
    0
    sage: try:
    ....:     exact_product(0^(k - 1), k, 0, 3)
    ....: except ValueError as err:
    ....:     print('ValueError: %s' % err)
    ValueError: the product 0^(k - 1) is not defined in the range

Parameter-linear factors over a concrete range are evaluated
directly, so that removable gamma poles do not become evaluation
errors::

    sage: exact_product(x + k, k, 0, 1)
    (x + 1)*x
    sage: exact_product(x + k, k, 0, 1).subs(x = -2)
    2

A vanishing product is an admissible value over a concrete range,
provided every factor is defined::

    sage: exact_product(k - 1, k, 0, 2)
    0
    sage: exact_product(0, k, 0, 3)
    0
    sage: exact_product((k - 1)*cos(k), k, 0, 2)
    0
    sage: try:
    ....:     exact_product(1/(k - 5), k, 0, 10)
    ....: except ValueError as err:
    ....:     print('ValueError: %s' % err)
    ValueError: the product 1/(k - 5) is not defined in the range

For a symbolic number of factors the gamma quotient is returned,
which then evaluates like any other gamma quotient, in particular
failing at the removable poles::

    sage: exact_product(x + k, k, 0, n)
    gamma(n + x + 1)/gamma(x)
    sage: try:
    ....:     exact_product(x + k, k, 0, n).subs(x = -2, n = 1)
    ....: except RuntimeError as err:
    ....:     print('RuntimeError: %s' % err)
    RuntimeError: indeterminate expression: 0 * infinity encountered.

The empty product is 1, whatever the body is::

    sage: exact_product(cos(k), k, 0, -1)
    1

The number of factors must be a nonnegative integer, or an
integer-linear expression in a single variable::

    sage: exact_product(x + k, k, 0, 1/2)
    Traceback (most recent call last):
    ...
    ValueError: the number of factors 3/2 is not a nonnegative integer
    sage: exact_product(x, k, 0, n^2)
    Traceback (most recent call last):
    ...
    ValueError: the number of factors n^2 + 1 is not linear in a single variable with rational coefficients

The residual part over a concrete range is evaluated directly::

    sage: exact_product((k + 1)*cos(k), k, 0, 1)
    2*cos(1)

The residual part is returned as an unevaluated product::

    sage: exact_product((k + 1)/(k + 2)*cos(k), k, 0, n)
    product(cos(k), k, 0, n)/(n + 2)

Linear factors with symbolic constant term are evaluated in terms of
gamma functions::

    sage: exact_product((x + k)/(k + 1), k, 0, n)
    gamma(n + x + 1)/(factorial(n + 1)*gamma(x))

Nothing can be done for genuinely non-rational bodies::

    sage: exact_product(cos(k), k, 0, n)
    Traceback (most recent call last):
    ...
    ValueError: the product cannot be evaluated

>>> from sage.all import *
>>> from sage.symbolic.gamma_products import exact_product
>>> k, n = var('k n')
>>> exact_product(Integer(2)**k, k, Integer(0), n)
2^(1/2*n^2 + 1/2*n)
>>> exact_product(factorial(k)/factorial(k + Integer(2)), k, Integer(0), n)
1/((n + 2)*factorial(n + 1)^2)
>>> exact_product((Integer(3)*k + Integer(1))*(Integer(3)*k + Integer(2))*factorial(k), k, Integer(0), n)
(1/3)^n*(3*n + 2)*(3*n + 1)*factorial(3*n)*product(factorial(k), k, 0, n)/factorial(n)

Powers with symbolic base and polynomial exponent, as well as
factors independent of the product variable, are evaluated; a zero
base is not aggregated, but evaluated at every position of a
concrete range, since it may have factors that are not defined::

>>> x = var('x')
>>> exact_product(x**k*(k + Integer(1)), k, Integer(0), n)
    x^(1/2*n^2 + 1/2*n)*factorial(n + 1)
>>> exact_product(sqrt(Integer(2))**k, k, Integer(0), n)
    2^(1/4*n^2 + 1/4*n)
>>> exact_product(Integer(0)**k, k, Integer(0), Integer(3))
    0
>>> try:
...     exact_product(Integer(0)**(k - Integer(1)), k, Integer(0), Integer(3))
... except ValueError as err:
...     print('ValueError: %s' % err)
    ValueError: the product 0^(k - 1) is not defined in the range

Parameter-linear factors over a concrete range are evaluated
directly, so that removable gamma poles do not become evaluation
errors::

>>> exact_product(x + k, k, Integer(0), Integer(1))
    (x + 1)*x
>>> exact_product(x + k, k, Integer(0), Integer(1)).subs(x = -Integer(2))
    2

A vanishing product is an admissible value over a concrete range,
provided every factor is defined::

>>> exact_product(k - Integer(1), k, Integer(0), Integer(2))
    0
>>> exact_product(Integer(0), k, Integer(0), Integer(3))
    0
>>> exact_product((k - Integer(1))*cos(k), k, Integer(0), Integer(2))
    0
>>> try:
...     exact_product(Integer(1)/(k - Integer(5)), k, Integer(0), Integer(10))
... except ValueError as err:
...     print('ValueError: %s' % err)
    ValueError: the product 1/(k - 5) is not defined in the range

For a symbolic number of factors the gamma quotient is returned,
which then evaluates like any other gamma quotient, in particular
failing at the removable poles::

>>> exact_product(x + k, k, Integer(0), n)
    gamma(n + x + 1)/gamma(x)
>>> try:
...     exact_product(x + k, k, Integer(0), n).subs(x = -Integer(2), n = Integer(1))
... except RuntimeError as err:
...     print('RuntimeError: %s' % err)
    RuntimeError: indeterminate expression: 0 * infinity encountered.

The empty product is 1, whatever the body is::

>>> exact_product(cos(k), k, Integer(0), -Integer(1))
    1

The number of factors must be a nonnegative integer, or an
integer-linear expression in a single variable::

>>> exact_product(x + k, k, Integer(0), Integer(1)/Integer(2))
    Traceback (most recent call last):
    ...
    ValueError: the number of factors 3/2 is not a nonnegative integer
>>> exact_product(x, k, Integer(0), n**Integer(2))
    Traceback (most recent call last):
    ...
    ValueError: the number of factors n^2 + 1 is not linear in a single variable with rational coefficients

The residual part over a concrete range is evaluated directly::

>>> exact_product((k + Integer(1))*cos(k), k, Integer(0), Integer(1))
    2*cos(1)

The residual part is returned as an unevaluated product::

>>> exact_product((k + Integer(1))/(k + Integer(2))*cos(k), k, Integer(0), n)
    product(cos(k), k, 0, n)/(n + 2)

Linear factors with symbolic constant term are evaluated in terms of
gamma functions::

>>> exact_product((x + k)/(k + Integer(1)), k, Integer(0), n)
    gamma(n + x + 1)/(factorial(n + 1)*gamma(x))

Nothing can be done for genuinely non-rational bodies::

>>> exact_product(cos(k), k, Integer(0), n)
    Traceback (most recent call last):
    ...
    ValueError: the product cannot be evaluated

from sage.symbolic.gamma_products import exact_product
k, n = var('k n')
exact_product(2^k, k, 0, n)
exact_product(factorial(k)/factorial(k + 2), k, 0, n)
exact_product((3*k + 1)*(3*k + 2)*factorial(k), k, 0, n)
x = var('x')
exact_product(x^k*(k + 1), k, 0, n)
exact_product(sqrt(2)^k, k, 0, n)
exact_product(0^k, k, 0, 3)
try:
    exact_product(0^(k - 1), k, 0, 3)
except ValueError as err:
    print('ValueError: %s' % err)
exact_product(x + k, k, 0, 1)
exact_product(x + k, k, 0, 1).subs(x = -2)
exact_product(k - 1, k, 0, 2)
exact_product(0, k, 0, 3)
exact_product((k - 1)*cos(k), k, 0, 2)
try:
    exact_product(1/(k - 5), k, 0, 10)
except ValueError as err:
    print('ValueError: %s' % err)
exact_product(x + k, k, 0, n)
try:
    exact_product(x + k, k, 0, n).subs(x = -2, n = 1)
except RuntimeError as err:
    print('RuntimeError: %s' % err)
exact_product(cos(k), k, 0, -1)
exact_product(x + k, k, 0, 1/2)
exact_product(x, k, 0, n^2)
exact_product((k + 1)*cos(k), k, 0, 1)
exact_product((k + 1)/(k + 2)*cos(k), k, 0, n)
exact_product((x + k)/(k + 1), k, 0, n)
exact_product(cos(k), k, 0, n)
sage.symbolic.gamma_products.factorial_product(f, var, a, b)[source]

Evaluate a finite product with a symbolic bound in factorial normal form.

INPUT:

  • f – a rational function of var with rational coefficients

  • var – the product variable

  • a – the lower bound, a rational number

  • b – the upper bound, such that the number of factors b - a + 1 is a nonnegative integer, or an integer-linear expression q*var + r in a single variable with \(q\) a positive integer and \(r\) a nonnegative integer

OUTPUT: the product \(\prod_{v=a}^{b} f(v)\), in the normal form described in simplify_gamma().

A ValueError is raised if the input is not of the required form.

EXAMPLES:

sage: from sage.symbolic.gamma_products import factorial_product
sage: k, n = var('k n')
sage: factorial_product(k, k, 1, n)
factorial(n)
sage: factorial_product(2*k + 1, k, 0, n - 1)
(1/2)^n*factorial(2*n)/factorial(n)
sage: factorial_product(1/(k + 1), k, 1, n)
1/factorial(n + 1)
sage: factorial_product(k^2, k, 1, 2*n)
factorial(2*n)^2

sage: f = (54*k^3 + 81*k^2 + 39*k + 6)/(2*k^3 + 9*k^2 + 13*k + 6)
sage: factorial_product(f, k, 0, n - 1)
factorial(3*n)/((2*n + 1)*(n + 1)*factorial(n)^3)

>>> from sage.all import *
>>> from sage.symbolic.gamma_products import factorial_product
>>> k, n = var('k n')
>>> factorial_product(k, k, Integer(1), n)
factorial(n)
>>> factorial_product(Integer(2)*k + Integer(1), k, Integer(0), n - Integer(1))
(1/2)^n*factorial(2*n)/factorial(n)
>>> factorial_product(Integer(1)/(k + Integer(1)), k, Integer(1), n)
1/factorial(n + 1)
>>> factorial_product(k**Integer(2), k, Integer(1), Integer(2)*n)
factorial(2*n)^2

>>> f = (Integer(54)*k**Integer(3) + Integer(81)*k**Integer(2) + Integer(39)*k + Integer(6))/(Integer(2)*k**Integer(3) + Integer(9)*k**Integer(2) + Integer(13)*k + Integer(6))
>>> factorial_product(f, k, Integer(0), n - Integer(1))
factorial(3*n)/((2*n + 1)*(n + 1)*factorial(n)^3)

from sage.symbolic.gamma_products import factorial_product
k, n = var('k n')
factorial_product(k, k, 1, n)
factorial_product(2*k + 1, k, 0, n - 1)
factorial_product(1/(k + 1), k, 1, n)
factorial_product(k^2, k, 1, 2*n)
f = (54*k^3 + 81*k^2 + 39*k + 6)/(2*k^3 + 9*k^2 + 13*k + 6)
factorial_product(f, k, 0, n - 1)

Finite products are evaluated exactly:

sage: factorial_product(k^2 - 1, k, 2, 6)
302400

>>> from sage.all import *
>>> factorial_product(k**Integer(2) - Integer(1), k, Integer(2), Integer(6))
302400

factorial_product(k^2 - 1, k, 2, 6)
sage.symbolic.gamma_products.simplify_gamma(expression, var)[source]

Rewrite a hypergeometric term as a product of factorials, in normal form.

INPUT:

  • expression – a symbolic expression which is a finite product and quotient of powers of rationals, of exponentials with rational base and rational-linear exponent, and of factorials, binomials, Pochhammer symbols and gamma functions with rational-linear argument. It must be a hypergeometric term in var and must be defined and nonzero at var = 0.

  • var – the variable

OUTPUT: an expression equivalent to expression of the form

\[C d^n \prod_i (A_i n + B_i)!^{e_i} \prod_k p_k(n)^{\varepsilon_k} \prod_j \Gamma(a_j n + b_j)^{\delta_j}\]

where \(C\) is a constant, \(d\) a nonzero rational, whose sign absorbs the sign of negative linear factors, \(A_i, B_i\) are nonnegative integers, \(p_k(n)\) are products of consecutive linear factors, and the residual gamma factors, if there are any, are in irreducible position: all shifts lie in the interval \([0, 1)\), shifts with the same slope \(a_j\) are distinct modulo the integers, and no application of the Gauss multiplication formula reduces the product further.

A ValueError is raised if expression is not a hypergeometric term, does not split into linear factors over QQ, or is zero or undefined at var = 0.

The normal form is not canonical: in particular, constants may be represented by gamma values at rational points or by evaluated powers of pi, and simplifying an already simplified expression again may change the representation.

EXAMPLES:

sage: from sage.symbolic.gamma_products import simplify_gamma
sage: n = var('n')
sage: simplify_gamma(binomial(3*n, n), n)
factorial(3*n)/(factorial(2*n)*factorial(n))
sage: simplify_gamma(gamma(n + 1/2), n)
sqrt(pi)*(1/4)^n*factorial(2*n)/factorial(n)
sage: simplify_gamma(factorial(n)*gamma(n + 1/2)/sqrt(pi), n)
(1/4)^n*factorial(2*n)

>>> from sage.all import *
>>> from sage.symbolic.gamma_products import simplify_gamma
>>> n = var('n')
>>> simplify_gamma(binomial(Integer(3)*n, n), n)
factorial(3*n)/(factorial(2*n)*factorial(n))
>>> simplify_gamma(gamma(n + Integer(1)/Integer(2)), n)
sqrt(pi)*(1/4)^n*factorial(2*n)/factorial(n)
>>> simplify_gamma(factorial(n)*gamma(n + Integer(1)/Integer(2))/sqrt(pi), n)
(1/4)^n*factorial(2*n)

from sage.symbolic.gamma_products import simplify_gamma
n = var('n')
simplify_gamma(binomial(3*n, n), n)
simplify_gamma(gamma(n + 1/2), n)
simplify_gamma(factorial(n)*gamma(n + 1/2)/sqrt(pi), n)

Purely exponential terms are accepted, in particular those whose ratio is constant:

sage: simplify_gamma(2^n, n)
2^n
sage: simplify_gamma((2/3)^n, n)
(2/3)^n

>>> from sage.all import *
>>> simplify_gamma(Integer(2)**n, n)
2^n
>>> simplify_gamma((Integer(2)/Integer(3))**n, n)
(2/3)^n

simplify_gamma(2^n, n)
simplify_gamma((2/3)^n, n)

Exponential atoms which are not factors of the term are kept, so that the conversion of the ratio to a rational function fails:

sage: simplify_gamma(2^n + 1, n)
Traceback (most recent call last):
...
ValueError: expression is not hypergeometric (ratio (2*2^n + 1)/(2^n + 1))

>>> from sage.all import *
>>> simplify_gamma(Integer(2)**n + Integer(1), n)
Traceback (most recent call last):
...
ValueError: expression is not hypergeometric (ratio (2*2^n + 1)/(2^n + 1))

simplify_gamma(2^n + 1, n)

A lone gamma with fractional shift cannot be reduced:

sage: simplify_gamma(gamma(n + 2/3)/gamma(2/3), n)
gamma(n + 2/3)/gamma(2/3)

>>> from sage.all import *
>>> simplify_gamma(gamma(n + Integer(2)/Integer(3))/gamma(Integer(2)/Integer(3)), n)
gamma(n + 2/3)/gamma(2/3)

simplify_gamma(gamma(n + 2/3)/gamma(2/3), n)

but two of them combine, by the Gauss multiplication formula:

sage: simplify_gamma(gamma(n + 1/3)*gamma(n + 2/3), n)
2/3*sqrt(3)*pi*(1/27)^n*factorial(3*n)/factorial(n)

When a complete reduction is impossible, the progress is committed
and the residual gamma factors are in irreducible position; here the
constants gamma(1/6) and gamma(5/6) admit no multiplication formula::

    sage: simplify_gamma(gamma(2*n + 1/3)*gamma(2*n + 2/3), n)
    4/3*sqrt(3)*pi^2*(1/729)^n*factorial(6*n)/(factorial(2*n)*gamma(5/6)*gamma(1/6))

>>> from sage.all import *
>>> simplify_gamma(gamma(n + Integer(1)/Integer(3))*gamma(n + Integer(2)/Integer(3)), n)
2/3*sqrt(3)*pi*(1/27)^n*factorial(3*n)/factorial(n)

When a complete reduction is impossible, the progress is committed
and the residual gamma factors are in irreducible position; here the
constants gamma(1/6) and gamma(5/6) admit no multiplication formula::

>>> simplify_gamma(gamma(Integer(2)*n + Integer(1)/Integer(3))*gamma(Integer(2)*n + Integer(2)/Integer(3)), n)
    4/3*sqrt(3)*pi^2*(1/729)^n*factorial(6*n)/(factorial(2*n)*gamma(5/6)*gamma(1/6))

simplify_gamma(gamma(n + 1/3)*gamma(n + 2/3), n)
simplify_gamma(gamma(2*n + 1/3)*gamma(2*n + 2/3), n)

The following term arises in the enumeration of certain trees; the denominator stems from a conversion of a product and was left in an unexpanded form:

sage: E = (3^(3*n)*gamma(n + 2/3)*gamma(n + 1/3)
....:      /((2*n*gamma(2/3)*gamma(1/3) + gamma(2/3)*gamma(1/3))
....:        *factorial(n + 1)*factorial(n)))
sage: simplify_gamma(E, n)
factorial(3*n)/((2*n + 1)*(n + 1)*factorial(n)^3)

>>> from sage.all import *
>>> E = (Integer(3)**(Integer(3)*n)*gamma(n + Integer(2)/Integer(3))*gamma(n + Integer(1)/Integer(3))
...      /((Integer(2)*n*gamma(Integer(2)/Integer(3))*gamma(Integer(1)/Integer(3)) + gamma(Integer(2)/Integer(3))*gamma(Integer(1)/Integer(3)))
...        *factorial(n + Integer(1))*factorial(n)))
>>> simplify_gamma(E, n)
factorial(3*n)/((2*n + 1)*(n + 1)*factorial(n)^3)

E = (3^(3*n)*gamma(n + 2/3)*gamma(n + 1/3)
     /((2*n*gamma(2/3)*gamma(1/3) + gamma(2/3)*gamma(1/3))
       *factorial(n + 1)*factorial(n)))
simplify_gamma(E, n)