# 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 :( – Stephen Jul 7 '17 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? – Skidsdev Jul 7 '17 at 14:15
• Possible duplicate of Calculate Power Series Coefficients – trichoplax Jul 7 '17 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 '17 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. – Peter Taylor Jul 8 '17 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++] & – ZaMoC Jul 7 '17 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. – Keyu Gan Jul 7 '17 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 InternalFromCoefficientList – ngenisis Jul 7 '17 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?) – Julian Wolf Jul 7 '17 at 18:51
• @ngenisis: What version are you using? On v10.0, your solution gives me Iterator {i,∞} does not have appropriate bounds. – Julian Wolf Jul 7 '17 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. – ZaMoC Jul 7 '17 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. – Julian Wolf Jul 7 '17 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. – Julian Wolf Jul 7 '17 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!