-- polynomial utils
a `nmul` n = map (*n) a
a `ndiv` n = map (`div` n) a

instance (Integral a) => Num [a] where
  (+) = zipWith (+)
  negate = map negate
  a * b = foldr f undefined b where
    f x z = (a `nmul` x) + (0 : z)
  abs _ = undefined
  signum _ = undefined
  fromInteger n = fromInteger n : repeat 0

-- replace x in polynomial with x^n
repl a n = concatMap (: replicate (n-1) 0) a

-- S2: (a^2 + b)/2
cycleIndexS2 a b = (a*a + b)`ndiv` 2

-- S4: (a^4 + 6 a^2 b + 8 a c + 3 b^2 + 6 d) / 24
cycleIndexS4 a b c d =	((a ^ 4) +
			 (a ^ 2 * b) `nmul` 6 +
			 (a * c) `nmul` 8 +
			 (b ^ 2) `nmul` 3 +
			 d `nmul` 6) `ndiv` 24


a598 = x1
-- A000598: A(x) = 1 + (1/6)*x*(A(x)^3 + 3*A(x)*A(x^2) + 2*A(x^3))
x1 = 1 : ((x1^3) + ((x2*x1)`nmul` 3) + (x3`nmul`2)) `ndiv` 6
x2 = x1`repl`2
x3 = x1`repl`3
x4 = x1`repl`4

-- A000678 = x CycleIndex(S4, A000598(x))
a678 = 0 : cycleIndexS4 x1 x2 x3 x4

-- A000599 = CycleIndex(S2, A000598(x) - 1)
a599 = cycleIndexS2 (0 : tail x1) (0 : tail x2)

-- A000602 = A000678(x) - A000599(x) + A000599(x^2)
a602 = a678 - a599 + x2

main = mapM_ print $ take 200 $ zip [0 ..] a602