The fractal plot modulo two Mathematica code. Clear[a, b, q, x, n, m, p, f]; p[x_, n_] := (x - 1)^(n + 1)* Sum[((-1)^(n + 1)*k^n)*x^k, {k, 0, Infinity}]/x; q[x_, n_] = D[p[x, n], {x, 2}]; f[n_] := CoefficientList[FullSimplify[ExpandAll[q[x, n]]], x]; a = Table[(f[n] + Reverse[f[n]])/4, {n, 3, 32}]; b = Table[If[m <= n, Mod[a[[n]][[m]], 2], 0], {m, 1, Length[a]}, {n, 1, Length[a]}]; ListDensityPlot[b, Mesh -> False]