Concrete Mathematics

Hypergeometric terms

The center stage, in recurrence solving and summations, play hypergeometric terms. Formally these are sequences annihilated by first order linear recurrence operators. In simple words if we are given term \(a(n)\) then it is hypergeometric if its consecutive term ratio is a rational function in \(n\).

To check if a sequence is of this type you can use the is_hypergeometric method which is available in Basic class. Here is simple example involving a polynomial:

>>> (n**2 + 1).is_hypergeometric(n)
True

Of course polynomials are hypergeometric but are there any more complicated sequences of this type? Here are some trivial examples:

>>> factorial(n).is_hypergeometric(n)
True
>>> binomial(n, k).is_hypergeometric(n)
True
>>> rf(n, k).is_hypergeometric(n)
True
>>> ff(n, k).is_hypergeometric(n)
True
>>> gamma(n).is_hypergeometric(n)
True
>>> (2**n).is_hypergeometric(n)
True

We see that all species used in summations and other parts of concrete mathematics are hypergeometric. Note also that binomial coefficients and both rising and falling factorials are hypergeometric in both their arguments:

>>> binomial(n, k).is_hypergeometric(k)
True
>>> rf(n, k).is_hypergeometric(k)
True
>>> ff(n, k).is_hypergeometric(k)
True

To say more, all previously shown examples are valid for integer linear arguments:

>>> factorial(2*n).is_hypergeometric(n)
True
>>> binomial(3*n+1, k).is_hypergeometric(n)
True
>>> rf(n+1, k-1).is_hypergeometric(n)
True
>>> ff(n-1, k+1).is_hypergeometric(n)
True
>>> gamma(5*n).is_hypergeometric(n)
True
>>> (2**(n-7)).is_hypergeometric(n)
True

However nonlinear arguments make those sequences fail to be hypergeometric:

>>> factorial(n**2).is_hypergeometric(n)
False
>>> (2**(n**3 + 1)).is_hypergeometric(n)
False

If not only the knowledge of being hypergeometric or not is needed, you can use hypersimp() function. It will try to simplify combinatorial expression and if the term given is hypergeometric it will return a quotient of polynomials of minimal degree. Otherwise is will return \(None\) to say that sequence is not hypergeometric:

>>> hypersimp(factorial(2*n), n)
2*(n + 1)*(2*n + 1)
>>> hypersimp(factorial(n**2), n)

Concrete Class Reference

class diofant.concrete.summations.Sum(function, *symbols, **assumptions)[source]

Represents an unevaluated summation.

Sum represents a finite or infinite series, with the first argument being the general form of terms in the series (which usually depend on the bound variable symbol), and the second argument being (symbol, start, end), with symbol taking all integer values from start through end (inclusive).

For start < end we adopt definition:

Sum(f(i), (i, start, end)) = -Sum(f(i), (i, end+1, start-1))

Notes

The summation convention for start < end described by Karr in [Kar81], especially definition 3 of section 1.4. The only difference with the reference is that Karr defines all sums with the upper limit being exclusive. This is in contrast to the usual mathematical notation, which we adopt, but does not affect the summation convention. Indeed we have:

\[\sum_{m \leq i < n} f_i = \sum_{i = m}^{n - 1} f_i\]

This convention allows us to preserve the splitting identity

\[\sum\limits_{i=m}^n f_i = \sum\limits_{i=m}^l f_i + \sum\limits_{i=l+1}^n f_i\]

regardless of the ordering of \(m\), \(l\) and \(n\).

Note that it also follows:

\[\sum\limits_{i=m}^{m-1} f_i = 0\]

Examples

>>> Sum(k**2, (k, 1, m)).doit()
m**3/3 + m**2/2 + m/6
>>> Sum(x**k/factorial(k), (k, 0, oo)).doit()
E**x

An example showing that the symbolic result of a summation is still valid for seemingly nonsensical values of the limits. Then the Karr convention allows us to give a perfectly valid interpretation to those sums by interchanging the limits according to the adopted rule:

>>> Sum(k, (k, 1, n)).doit()
n**2/2 + n/2
>>> _.subs({n: -4})
6
>>> Sum(-n, (n, -3, 0)).doit()
6

References

euler_maclaurin(m=0, n=0, eps=0, eval_integral=True)[source]

Return an Euler-Maclaurin approximation of self, where m is the number of leading terms to sum directly and n is the number of terms in the tail.

With m = n = 0, this is simply the corresponding integral plus a first-order endpoint correction.

Returns (s, e) where s is the Euler-Maclaurin approximation and e is the estimated error (taken to be the magnitude of the first omitted term in the tail):

>>> Sum(1/k, (k, 2, 5)).doit().evalf()
1.28333333333333
>>> s, e = Sum(1/k, (k, 2, 5)).euler_maclaurin()
>>> s
-log(2) + 7/20 + log(5)
>>> print(sstr((s.evalf(), e.evalf()), full_prec=True))
(1.26629073187415, 0.0175000000000000)

The endpoints may be symbolic:

>>> s, e = Sum(1/k, (k, a, b)).euler_maclaurin()
>>> s
-log(a) + log(b) + 1/(2*b) + 1/(2*a)
>>> e
Abs(1/(12*b**2) - 1/(12*a**2))

If the function is a polynomial of degree at most 2n+1, the Euler-Maclaurin formula becomes exact (and e = 0 is returned):

>>> Sum(k, (k, 2, b)).euler_maclaurin()
(b**2/2 + b/2 - 1, 0)
>>> Sum(k, (k, 2, b)).doit()
b**2/2 + b/2 - 1

With a nonzero \(eps\) specified, the summation is ended as soon as the remainder term is less than the epsilon.

findrecur(F=Function('F'), n=None)[source]

Find a recurrence formula for the summand of the sum.

Given a sum \(f(n) = \sum_k F(n, k)\), where \(F(n, k)\) is doubly hypergeometric (that’s, both \(F(n + 1, k)/F(n, k)\) and \(F(n, k + 1)/F(n, k)\) are rational functions of \(n\) and \(k\)), we find a recurrence for the summand \(F(n, k)\) of the form

\[\sum_{i=0}^I\sum_{j=0}^J a_{i,j}F(n - j, k - i) = 0\]

Examples

>>> s = Sum(factorial(n)/(factorial(k)*factorial(n - k)), (k, 0, oo))
>>> s.findrecur()
-F(n, k) + F(n - 1, k) + F(n - 1, k - 1)

Notes

We use Sister Celine’s algorithm, see [PetkovvsekWZ97], Ch. 4.

reverse_order(*indices)[source]

Reverse the order of a limit in a Sum.

Parameters:

*indices (list) – The selectors in the argument indices specify some indices whose limits get reversed. These selectors are either variable names or numerical indices counted starting from the inner-most limit tuple.

Examples

>>> Sum(x, (x, 0, 3)).reverse_order(x)
Sum(-x, (x, 4, -1))
>>> Sum(x*y, (x, 1, 5), (y, 0, 6)).reverse_order(x, y)
Sum(x*y, (x, 6, 0), (y, 7, -1))
>>> Sum(x, (x, a, b)).reverse_order(x)
Sum(-x, (x, b + 1, a - 1))
>>> Sum(x, (x, a, b)).reverse_order(0)
Sum(-x, (x, b + 1, a - 1))

While one should prefer variable names when specifying which limits to reverse, the index counting notation comes in handy in case there are several symbols with the same name.

>>> s = Sum(x**2, (x, a, b), (x, c, d))
>>> s
Sum(x**2, (x, a, b), (x, c, d))
>>> s0 = s.reverse_order(0)
>>> s0
Sum(-x**2, (x, b + 1, a - 1), (x, c, d))
>>> s1 = s0.reverse_order(1)
>>> s1
Sum(x**2, (x, b + 1, a - 1), (x, d + 1, c - 1))

Of course we can mix both notations:

>>> Sum(x*y, (x, a, b), (y, 2, 5)).reverse_order(x, 1)
Sum(x*y, (x, b + 1, a - 1), (y, 6, 1))
>>> Sum(x*y, (x, a, b), (y, 2, 5)).reverse_order(y, x)
Sum(x*y, (x, b + 1, a - 1), (y, 6, 1))

References

class diofant.concrete.products.Product(function, *symbols, **assumptions)[source]

Represents an unevaluated products.

Product represents a finite or infinite product, with the first argument being the general form of terms in the series (which usually depend on the bound variable symbol), and the second argument being (symbol, start, end), with symbol taking all integer values from start through end (inclusive).

Notes

We follow the the analogue of the summation convention described by Karr [Kar81], adopted by the Sum:

\[\prod\limits_{i=m}^n f_i = \frac{1}{\prod\limits_{i=n+1}^{m-1}f_i}\]

Examples

>>> Product(k**2, (k, 1, m)).doit()
factorial(m)**2

Products with the lower limit being larger than the upper one:

>>> Product(1/k, (k, 6, 1)).doit()
120
>>> Product(k, (k, 2, 5)).doit()
120

The empty product:

>>> Product(k, (k, n, n-1)).doit()
1

References

reverse_order(*indices)[source]

Reverse the order of a limit in a Product.

Parameters:

*indices (list) – The selectors in the argument indices specify some indices whose limits get reversed. These selectors are either variable names or numerical indices counted starting from the inner-most limit tuple.

Examples

>>> P = Product(x, (x, a, b))
>>> Pr = P.reverse_order(x)
>>> Pr
Product(1/x, (x, b + 1, a - 1))
>>> Pr = Pr.doit()
>>> Pr
1/RisingFactorial(b + 1, a - b - 1)
>>> simplify(Pr)
gamma(b + 1)/gamma(a)
>>> P = P.doit()
>>> P
RisingFactorial(a, -a + b + 1)
>>> simplify(P)
gamma(b + 1)/gamma(a)

While one should prefer variable names when specifying which limits to reverse, the index counting notation comes in handy in case there are several symbols with the same name.

>>> s = Sum(x*y, (x, a, b), (y, c, d))
>>> s
Sum(x*y, (x, a, b), (y, c, d))
>>> s0 = s.reverse_order(0)
>>> s0
Sum(-x*y, (x, b + 1, a - 1), (y, c, d))
>>> s1 = s0.reverse_order(1)
>>> s1
Sum(x*y, (x, b + 1, a - 1), (y, d + 1, c - 1))

Of course we can mix both notations:

>>> Sum(x*y, (x, a, b), (y, 2, 5)).reverse_order(x, 1)
Sum(x*y, (x, b + 1, a - 1), (y, 6, 1))
>>> Sum(x*y, (x, a, b), (y, 2, 5)).reverse_order(y, x)
Sum(x*y, (x, b + 1, a - 1), (y, 6, 1))

References

class diofant.concrete.expr_with_limits.ExprWithLimits(function, *symbols, **assumptions)[source]

Represents an expression with limits.

as_dummy()[source]

Replace instances of the given dummy variables with explicit dummy counterparts to make clear what are dummy variables and what are real-world symbols in an object.

Examples

>>> Integral(x, (x, x, y), (y, x, y)).as_dummy()
Integral(_x, (_x, x, _y), (_y, x, y))

If the object supports the “integral at” limit (x,) it is not treated as a dummy, but the explicit form, (x, x) of length 2 does treat the variable as a dummy.

>>> Integral(x, x).as_dummy()
Integral(x, x)
>>> Integral(x, (x, x)).as_dummy()
Integral(_x, (_x, x))

If there were no dummies in the original expression, then the the symbols which cannot be changed by subs() are clearly seen as those with an underscore prefix.

See also

diofant.concrete.expr_with_limits.ExprWithLimits.variables

Lists the integration variables

property free_symbols

This method returns the symbols in the object, excluding those that take on a specific value (i.e. the dummy symbols).

Examples

>>> Sum(x, (x, y, 1)).free_symbols
{y}
property function

Return the function applied across limits.

Examples

>>> Integral(x**2, x).function
x**2
property is_number

Return True if the Sum has no free symbols, else False.

property limits

Return the limits of expression.

Examples

>>> from diofant.abc import i
>>> Integral(x**i, (i, 1, 3)).limits
((i, 1, 3),)
property variables

Return a list of the dummy variables

>>> from diofant.abc import i
>>> Sum(x**i, (i, 1, 3)).variables
[i]
class diofant.concrete.expr_with_intlimits.ExprWithIntLimits(function, *symbols, **assumptions)[source]

Represents an expression with integer limits.

change_index(var, trafo, newvar=None)[source]

Change index of a Sum or Product.

Perform a linear transformation \(x \mapsto a x + b\) on the index variable \(x\). For \(a\) the only values allowed are \(\pm 1\). A new variable to be used after the change of index can also be specified.

Parameters:
  • var (Symbol) – specifies the index variable \(x\) to transform.

  • trafo (Expr) – The linear transformation in terms of var.

  • newvar (Symbol, optional) – Replacement symbol to be used instead of var in the final expression.

Examples

>>> from diofant.abc import i, j, l, u, v
>>> s = Sum(x, (x, a, b))
>>> s.doit()
-a**2/2 + a/2 + b**2/2 + b/2
>>> sn = s.change_index(x, x + 1, y)
>>> sn
Sum(y - 1, (y, a + 1, b + 1))
>>> sn.doit()
-a**2/2 + a/2 + b**2/2 + b/2
>>> sn = s.change_index(x, -x, y)
>>> sn
Sum(-y, (y, -b, -a))
>>> sn.doit()
-a**2/2 + a/2 + b**2/2 + b/2
>>> sn = s.change_index(x, x+u)
>>> sn
Sum(-u + x, (x, a + u, b + u))
>>> sn.doit()
-a**2/2 - a*u + a/2 + b**2/2 + b*u + b/2 - u*(-a + b + 1) + u
>>> simplify(sn.doit())
-a**2/2 + a/2 + b**2/2 + b/2
>>> sn = s.change_index(x, -x - u, y)
>>> sn
Sum(-u - y, (y, -b - u, -a - u))
>>> sn.doit()
-a**2/2 - a*u + a/2 + b**2/2 + b*u + b/2 - u*(-a + b + 1) + u
>>> simplify(sn.doit())
-a**2/2 + a/2 + b**2/2 + b/2
>>> p = Product(i*j**2, (i, a, b), (j, c, d))
>>> p
Product(i*j**2, (i, a, b), (j, c, d))
>>> p2 = p.change_index(i, i+3, k)
>>> p2
Product(j**2*(k - 3), (k, a + 3, b + 3), (j, c, d))
>>> p3 = p2.change_index(j, -j, l)
>>> p3
Product(l**2*(k - 3), (k, a + 3, b + 3), (l, -d, -c))

When dealing with symbols only, we can make a general linear transformation:

>>> sn = s.change_index(x, u*x+v, y)
>>> sn
Sum((-v + y)/u, (y, b*u + v, a*u + v))
>>> sn.doit()
-v*(a*u - b*u + 1)/u + (a**2*u**2/2 + a*u*v + a*u/2 - b**2*u**2/2 - b*u*v + b*u/2 + v)/u
>>> simplify(sn.doit())
a**2*u/2 + a/2 - b**2*u/2 + b/2

However, the last result can be inconsistent with usual summation where the index increment is always 1. This is obvious as we get back the original value only for u equal +1 or -1.

index(x)[source]

Return the index of a dummy variable in the list of limits.

Note that we start counting with 0 at the inner-most limits tuple.

Parameters:

x (Symbol) – a dummy variable

Examples

>>> Sum(x*y, (x, a, b), (y, c, d)).index(x)
0
>>> Sum(x*y, (x, a, b), (y, c, d)).index(y)
1
>>> Product(x*y, (x, a, b), (y, c, d)).index(x)
0
>>> Product(x*y, (x, a, b), (y, c, d)).index(y)
1
reorder(*arg)[source]

Reorder limits in a expression containing a Sum or a Product.

Parameters:

*arg (list of tuples) – These tuples can contain numerical indices or index variable names or involve both.

Examples

>>> from diofant.abc import e, f
>>> Sum(x*y, (x, a, b), (y, c, d)).reorder((x, y))
Sum(x*y, (y, c, d), (x, a, b))
>>> Sum(x*y*z, (x, a, b), (y, c, d), (z, e, f)).reorder((x, y), (x, z), (y, z))
Sum(x*y*z, (z, e, f), (y, c, d), (x, a, b))
>>> P = Product(x*y*z, (x, a, b), (y, c, d), (z, e, f))
>>> P.reorder((x, y), (x, z), (y, z))
Product(x*y*z, (z, e, f), (y, c, d), (x, a, b))

We can also select the index variables by counting them, starting with the inner-most one:

>>> Sum(x**2, (x, a, b), (x, c, d)).reorder((0, 1))
Sum(x**2, (x, c, d), (x, a, b))

And of course we can mix both schemes:

>>> Sum(x*y, (x, a, b), (y, c, d)).reorder((y, x))
Sum(x*y, (y, c, d), (x, a, b))
>>> Sum(x*y, (x, a, b), (y, c, d)).reorder((y, 0))
Sum(x*y, (y, c, d), (x, a, b))
reorder_limit(x, y)[source]

Interchange two limit tuples of a Sum or Product expression.

Parameters:

x, y (int) – are integers corresponding to the index variables of the two limits which are to be interchanged.

Examples

>>> from diofant.abc import e, f
>>> Sum(x*y*z, (x, a, b), (y, c, d), (z, e, f)).reorder_limit(0, 2)
Sum(x*y*z, (z, e, f), (y, c, d), (x, a, b))
>>> Sum(x**2, (x, a, b), (x, c, d)).reorder_limit(1, 0)
Sum(x**2, (x, c, d), (x, a, b))
>>> Product(x*y*z, (x, a, b), (y, c, d), (z, e, f)).reorder_limit(0, 2)
Product(x*y*z, (z, e, f), (y, c, d), (x, a, b))

Concrete Functions Reference

diofant.concrete.summations.summation(f, *symbols, **kwargs)[source]

Compute the summation of f with respect to symbols.

The notation for symbols is similar to the notation used in Integral. summation(f, (i, a, b)) computes the sum of f with respect to i from a to b, i.e.,

                            b
                          ____
                          \   `
summation(f, (i, a, b)) =  )    f
                          /___,
                          i = a

If it cannot compute the sum, it returns an unevaluated Sum object. Repeated sums can be computed by introducing additional symbols tuples:

>>> i = symbols('i', integer=True)
>>> summation(2*i - 1, (i, 1, n))
n**2
>>> summation(1/2**i, (i, 0, oo))
2
>>> summation(1/log(n)**n, (n, 2, oo))
Sum(log(n)**(-n), (n, 2, oo))
>>> summation(i, (i, 0, n), (n, 0, m))
m**3/6 + m**2/2 + m/3
>>> summation(x**n/factorial(n), (n, 0, oo))
E**x
diofant.concrete.products.product(*args, **kwargs)[source]

Compute the product.

The notation for symbols is similar to the notation used in Sum or Integral. product(f, (i, a, b)) computes the product of f with respect to i from a to b, i.e.,

                             b
                           _____
product(f(n), (i, a, b)) = |   | f(n)
                           |   |
                           i = a

If it cannot compute the product, it returns an unevaluated Product object. Repeated products can be computed by introducing additional symbols tuples:

>>> i = symbols('i', integer=True)
>>> product(i, (i, 1, k))
factorial(k)
>>> product(m, (i, 1, k))
m**k
>>> product(i, (i, 1, k), (k, 1, n))
Product(factorial(k), (k, 1, n))
diofant.concrete.gosper.gosper_normal(f, g, n)[source]

Compute the Gosper’s normal form of f and g.

Given relatively prime univariate polynomials f and g, rewrite their quotient to a normal form defined as follows:

\[\frac{f(n)}{g(n)} = Z \cdot \frac{A(n) C(n+1)}{B(n) C(n)}\]

where Z is an arbitrary constant and A, B, C are monic polynomials in n with the following properties:

  1. \(\gcd(A(n), B(n+h)) = 1 \forall h \in \mathbb{N}\)

  2. \(\gcd(B(n), C(n+1)) = 1\)

  3. \(\gcd(A(n), C(n)) = 1\)

This normal form, or rational factorization in other words, is a crucial step in Gosper’s algorithm and in solving of difference equations. It can be also used to decide if two hypergeometric terms are similar or not.

This procedure will return a tuple containing elements of this factorization in the form (Z*A, B, C).

Examples

>>> gosper_normal(4*n + 5, 2*(4*n + 1)*(2*n + 3), n)
(Poly(1/4, n, domain='QQ'), Poly(n + 3/2, n, domain='QQ'),
 Poly(n + 1/4, n, domain='QQ'))
diofant.concrete.gosper.gosper_term(f, n)[source]

Compute Gosper’s hypergeometric term for f.

Suppose f is a hypergeometric term such that:

\[s_n = \sum_{k=0}^{n-1} f_k\]

and \(f_k\) doesn’t depend on \(n\). Returns a hypergeometric term \(g_n\) such that \(g_{n+1} - g_n = f_n\).

Examples

>>> gosper_term((4*n + 1)*factorial(n)/factorial(2*n + 1), n)
(-n - 1/2)/(n + 1/4)
diofant.concrete.gosper.gosper_sum(f, k)[source]

Gosper’s hypergeometric summation algorithm.

Given a hypergeometric term f such that:

\[s_n = \sum_{k=0}^{n-1} f_k\]

and \(f(n)\) doesn’t depend on \(n\), returns \(g_{n} - g(0)\) where \(g_{n+1} - g_n = f_n\), or None if \(s_n\) can not be expressed in closed form as a sum of hypergeometric terms.

Examples

>>> gosper_sum((4*k + 1)*factorial(k)/factorial(2*k + 1), (k, 0, n))
(-factorial(n) + 2*factorial(2*n + 1))/factorial(2*n + 1)

References