# Find the coefficients of a rational generating function

If we write a sequence of numbers as the coefficients of a power series, then that power series is called the (ordinary) generating function (or G.f.) of that sequence. That is, if for some function F(x) and series of integers a(n) we have:

a(0) + a(1)x + a(2)x^2 + a(3)x^3 + a(4)x^4 + ... = F(x)

Then F(x) is the generating function of a. For example, the geometric series tells us that:

1 + x + x^2 + x^3 + x^4 + ... = 1/(1-x)

So the generating function of 1, 1, 1, ... is 1/(1-x). If we differentiate both sides of the equation above and multiply by x we get the following equality:

x + 2x^2 + 3x^3 + 4x^4 + ... = x/(1-x)^2

So the generating function of 1, 2, 3, ... is x/(1-x)^2. Generating functions are a very powerful tool, and you can do many useful things with them. A short introduction can be found here, but for a really thorough explanation there is the amazing book generatingfunctionology.

In this challenge you will take a rational function (the quotient of two polynomials with integer coefficients) as input as two arrays of integer coefficients, first the numerator then the denominator. For example the function f(x) = x / (1 - x - x^2) will be encoded as [0, 1], [1, -1, -1] in the input.

Given this input your program must infinitely print the coefficients of the power series that equals the generating function, one per line, starting at the coefficient of x, then x^2, etc.

Examples:

[1], [1, -1] -> 1, 1, 1, 1, 1, 1, 1, ...
[1], [2, -2] -> 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, ...
[0, 1], [1, -2, 1] -> 1, 2, 3, 4, 5, 6, 7, 8, ...
[0, 1], [1, -1, -1] -> 1, 1, 2, 3, 5, 8, 13, 21, 34, ...
[1], [1, -2] -> 1, 2, 4, 8, 16, 32, 64, 128, ...
[0, 1, 1], [1, -3, 3, -1] -> 1, 4, 9, 16, 25, 36, ...
• Crap, my language is built for sequences this, but I can't really do multidimensional array input :( Jul 7, 2017 at 14:02
• I'm really just not mathematically-minded enough for this spec, any chance you could post more of a layman's explanation for us common folk? Jul 7, 2017 at 14:15
• Possible duplicate of Calculate Power Series Coefficients Jul 7, 2017 at 23:39
• @trichoplax That always forces the numerator to be 1, which is not the same. For example it cannot express my last example, the squares.
– orlp
Jul 8, 2017 at 5:53
• An alternative way of phrasing this is that it evaluates a general linear recurence. In that way it generalises this question, and might serve as a dupe target for future recurrence questions. Jul 8, 2017 at 8:58

z=0:z
(a:b)%y@(c:d)=a/c:zipWith(-)(b++z)(map(a/c*)d++z)%y
_%_=z

Try it online!

Defines an operator % returning an infinite lazy list of coefficients. The list is zero-indexed, so the constant coefficient is included.

# Mathematica, 64 8390 bytes

Do[Echo@Limit[D[#/#2/i!&@@Fold[x#+#2&]/@#,{x,i}],x->0],{i,∞}‌​]&

Thanks to @ngenisis and @Jenny_mathy !

Take input as two lists.

Need Alt+. to terminate the execution to see the result. Frontend may crash due to rapid output.

83 bytes version (@Jenny_mathy):

i=1;v=Tr[#*x^Range@Length@#]&;While[1<2,Echo@Limit[D[v@#/v@#2/i!,{x,i}],x->0];i++]&
• 83 bytes: i = 1; v = Tr[#*x^Range@Length@#] &; While[1 < 2, Echo@Limit[D[v@#/v@#2/i!, {x, i}], x -> 0]; i++] & Jul 7, 2017 at 17:09
• @Jenny_mathy Sorry for bothering. I figure it out that there's some junk invisible Unicode characters in your first comment. Once cleaned up, the code is OK. Jul 7, 2017 at 17:15
• 64 bytes: Do[Echo@Limit[D[#/#2/i!&@@Fold[x#+#2&]/@#,{x,i}],x->0],{i,∞}]&. This assumes the input is a list of two lists and the coefficients are in order of descending degree. The only built-in I know of to do what v does is Internal`FromCoefficientList
– user61980
Jul 7, 2017 at 17:55
• Does running this repeatedly work? I think a couple of extra parentheses might be necessary to put i inside the lambda. (On the other hand, I'm not really sure whether the ability to run repeatedly is relevant when the goal is to print an infinite list…has there been a meta consensus on this?) Jul 7, 2017 at 18:51
• @ngenisis: What version are you using? On v10.0, your solution gives me Iterator {i,∞} does not have appropriate bounds. Jul 7, 2017 at 19:52

## CJam (22 bytes)

{\{(W$(@\/_pW*f*.+1}g} Online demo. Note that as many of the existing answers, this includes the 0th coefficient in the output. ### Dissection { e# Define a block which takes numerator N and denominator D as arguments \ e# Flip to put D at the bottom, since that won't change { e# Infinite loop: ( e# Pop the first element of (the modified) N W$(     e#   Copy D and pop its first element
e#   Stack: D N[1:] N[0] D[1:] D[0]
@\/     e#   Bring N[0] to top, flip, divide
e#   Stack: D N[1:] D[1:] N[0]/D[0]
_p      e#   Print a copy
W*f*.+  e#   Multiply by -1, multiply all, pointwise add
e#   Stack: D N[1:]-(N[0]/D[0])*D[1:]
1}g
}

# Mathematica, 86 79 bytes

f=x^Range@Length@#.#&;For[n=1,8>3,Print@SeriesCoefficient[f@#/f@#2,{x,0,n++}]]&

Takes input as two separate lists (numerator coefficients, denominator coefficients). If the input can be taken directly as a fraction of polynomials rather than as lists of coefficients, this can be shortened significantly.

It seems that Do can work with infinite bounds in v11. I can't test this locally, but, if this is the case, then this solution can be shortened to 75 bytes:

f=x^Range@Length@#.#&;Do[Print@SeriesCoefficient[f@#/f@#2,{x,0,n}],{n,∞}]&
• last test case doesn't begin with 0. Jul 7, 2017 at 18:56
• @Jenny_mathy: shoot, thanks for the heads up. Looks like the test cases expect starting from the first in stead of the zeroth…pretty sure this should let me save a few bytes. Jul 7, 2017 at 19:05
• @Jenny_mathy: I think the test cases might be wonky. Starting n from 1 in stead of 0, this gives the same results as your solution; both fail, though, on the second to last test case, which this solution passes when starting n from 0. Jul 7, 2017 at 20:00

# Pyth, 23 bytes

JE#
KchQhJ=t-M.t,Q*LKJ0

Try it online!

### How it works

Q = eval(input())
JE                     J = eval(input())
#                    infinite loop:
chQhJ                   Q[0]/J[0]
K                        assign that to K (and print it, because of the preceding newline)
*LKJ       K times every element of J
,Q           [Q, that]
.t      0      transpose, padding with 0s
-M               subtract each pair
t                 remove the first element
=                  assign that back to Q

# Pari/GP, 66 bytes

p->q->for(i=0,oo,print(Pol(Polrev(p)/(Polrev(q)+O(x*x^i)))\x^i%x))

Try it online!