The fractal plot modulo two Mathematica code.


Clear[a, b, q, x, n, m, p, f];
p[x_, n_]:= 2^n*(1-x)^(n+1)* LerchPhi[x, -n, 1/2];
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, 2, 32}]; 
b = Table[If[m <= n, Mod[a[[n]][[m]], 2], 0], {m, 1, Length[a]}, {n, 1, Length[a]}];
ListDensityPlot[b, Mesh -> False]