4

I have a rational function: f(x) = P(x)/Q(x). For example:

f(x) = (5x + 3)/(1-x^2)

Because f(x) is a generating function it can be written as:

f(x) = a0 + a1*x + a2*x² + ... + a_n*x^n + ... = P(x)/Q(x)

How can I use sympy to find the nth term of the generating function f(x) (that is a_n)?

If there is no such implementation in Sympy, I am curious also to know if this implemented in other packages, such as Maxima.

I appreciate any help.

JohanC
  • 71,591
  • 8
  • 33
  • 66
  • Did you look at https://stackoverflow.com/questions/46422538/how-to-solve-an-algebraic-equation-in-formal-power-series ? – JohanC Dec 19 '19 at 17:54

5 Answers5

6

To get the general formula for a_n of the generating function of a rational form , SymPy's rational_algorithm can be used. For example:

from sympy import simplify
from sympy.abc import x, n
from sympy.series.formal import rational_algorithm

f = (5*x + 3)/(1-x**2)
func_n, independent_term, order = rational_algorithm(f, x, n, full=True)
print(f"The general formula for a_n is {func_n}")
for k in range(10):
    print(f"a_{k} = {simplify(func_n.subs(n, k))}")

Output:

The general formula for a_n is (-1)**(-n - 1) + 4
a_0 = 3
a_1 = 5
a_2 = 3
a_3 = 5
a_4 = 3
a_5 = 5
a_6 = 3
a_7 = 5
a_8 = 3
a_9 = 5

Here is another example:

f = x / (1 - x - 2 * x ** 2)
func_n, independent_term, order = rational_algorithm(f, x, n, full=True)
print(f"The general formula for a_n is {func_n.simplify()}")
print("First terms:", [simplify(func_n.subs(n, k)) for k in range(20)])

The general formula for a_n is 2**n/3 - (-1)**(-n)/3
First terms: [0, 1, 1, 3, 5, 11, 21, 43, 85, 171, 341, 683, 1365, 2731, 5461, 10923, 21845, 43691, 87381, 174763]
JohanC
  • 71,591
  • 8
  • 33
  • 66
1

You could take the kth derivative and substitute 0 for x and divide by factorial(k):

>>> f = (5*x + 3) / (1-x**2)
>>> f.diff(x, 20).subs(x, 0)/factorial(20)
3

The reference here talks about rational generating functions. Looking for a recurrence you can see the pattern pretty quickly using differentiation:

[f.diff(x,i).subs(x,0)/factorial(i) for i in range(6)]
[3, 5, 3, 5, 3, 5]
smichr
  • 16,948
  • 2
  • 27
  • 34
0

Adapting the approach of this post, you could try the following:

from sympy import *
from sympy.abc import x

f = (5*x + 3) / (1-x**2)
print(f.series(n=20))
k = 50
coeff50 = Poly(series(f, x, n=k + 1).removeO(), x).coeff_monomial(x ** k)
print(f"The coeffcient of x^{k} of the generating function of {f} is {coeff50}")

# to get the first 100 coeffcients (reversing the list to get a[0] the
#   coefficient of x**0 etc.):
a = Poly(series(f, x, n=100).removeO(), x).all_coeffs()[::-1]

Output:

3 + 5*x + 3*x**2 + 5*x**3 + 3*x**4 + 5*x**5 + 3*x**6 + 5*x**7 + 3*x**8 + 5*x**9 + 3*x**10 + 5*x**11 + 3*x**12 + 5*x**13 + 3*x**14 + 5*x**15 + 3*x**16 + 5*x**17 + 3*x**18 + 5*x**19 + O(x**20)
The coeffcient of x^50 of the generating function of (5*x + 3)/(1 - x**2) is 3

Following this example at Cut The Knot, the approach can be used to find out the number of ways an amount n can be paid with coins of 1, 5, 10, 25 and 50 cents.

f = 1/((1 - x)*(1 - x**5)*(1 - x**10)*(1 - x**25)*(1 - x**50))
a = Poly(series(f, x, n=101).removeO(), x).all_coeffs()[::-1]
print(a[50]) # there are 50 ways to pay 50 cents
print(a[100]) # there are 292 ways to pay 100 cents
JohanC
  • 71,591
  • 8
  • 33
  • 66
0

In maxima:

powerseries((5*x+3)/(1-x^2),x,0);

returns

powerseries_maxima_result

Use part to extract the generator:

part(''%,1);
(4-(-1)^i1)x^i1

and coeff to get the coefficient:

a(i1) := coeff(''%, x, i1);
[a(0), a(1), a(2)];
[3, 5, 3]
Eelvex
  • 8,975
  • 26
  • 42
0

Another nice way to approach this is to use the ring series:

>>> from sympy.polys.ring_series import rs_mul, rs_pow
>>> from sympy.polys.rings import ring
>>> R,x=ring('x', ZZ)
>>> from sympy import ZZ
>>> R,x=ring('x', ZZ)
>>> nmax = 100
>>> s = rs_mul(5*x+3, rs_pow(1-x**2, -1, x, nmax+1), x, nmax+1)
>>> [s.coeff(x**i) for i in (2,3,5,17,100)]
[3, 5, 5, 5, 3]
smichr
  • 16,948
  • 2
  • 27
  • 34