Find the coefficients of a rational generating function

Question

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, ...

Answer 1

Haskell, 63 bytes

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.

Answer 2

Mathematica, 64 83 90 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++]&

Answer 3

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
}

Answer 4

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,∞}]&

Answer 5

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

Answer 6

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!