45 lines
1.1 KiB
Haskell
45 lines
1.1 KiB
Haskell
-- 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
|