原本是想直接续在这帖：https://tieba.baidu.com/p/2134509239的后面直接回复的，但是考虑到那帖楼层数已经较多，默默发在那边无法满足我的虚荣心（殴）所以决定重开一帖。（其实@隨意超 似乎已经在那帖11楼找对了方向，但却似乎没有引起后续的参与者的重视。）
先上代码和图再讲解吧。2GB内存2GHz老爷机，3 min出图：
modifiedroots[c_List] :=
Block[{a = DiagonalMatrix[ConstantArray[1., Length@c - 1], -1]},
a[[1]] = -c;
Eigenvalues[a]];
s = 1000;
l = ConstantArray[0., {s + 1, s + 1}];
l[[#, #2]] += 1. & @@@
Round[1 + s Rescale@
Function[z, {Im@z, Re@z}, Listable]@
Flatten[modifiedroots /@ Tuples[{-1., 1.}, 18]]]; // AbsoluteTiming
ArrayPlot[UnitStep[80 - l] l, ColorFunction -> "AvocadoColors"]
当然，我说“复现”其实是标题党了，因为原图算到了24次多项式，我这幅图只算到了19次。（这是2G内存的极限，次数每加1，内存占用大约要翻一番。）但是如大家所见，这幅图的效果已经和那幅Sam Derbyshire画了3天3夜的图相差无几了。
下面简单讲解一下这段代码。要做到在多数人可以
接受的时间内使用多数人
可以接受的内存获得
和那幅图相似的效果，需要：
1　省略对这个特别的问题而言其实不必要的（主要是符号）计算。
2　昨晚对像素点十分有限的屏幕而言其实不必要的仅有着细微差别的数据。
3　正确地注意到原图的颜色代表的是什么。
先说第一点。因为这幅图表示的是幂方程的根，所以我们自然需要求解幂方程。如何解方程？我们的第一反应是Solve。是的，Solve可以求解这个问题，但是，作为一个通用求解器，它对这个问题而言太慢了。有没有更为专用的函数可以求解这个问题？在本文开头的帖子里@草红样 找到了NRoots，一个专解多项式方程的函数，从而大大提速了求解。可是，这是否是终点？不，NRoots依旧需要进行多项式结构的分析，我们能不能把这一部也给省掉，直接使用（以确定次序排列的）多项式系数表作为输入来求根？再一次，答案是肯定的，那就是我们上面使用的modifiedroots。（对于NRoots的提速问题的讨论其实尚未结束，有兴趣的参考这帖：https://tieba.baidu.com/p/3575212088）
另，不难注意到“全体系数为正负1的全体n次幂方程”其实有一半是重复的。上面的modifiedroots已经利用了这个特性。
再说后两点。（这两点其实是相联系的，所以一起讲。）对方程的求解将会产生大量的数据，这倒是其次，真正够呛的是如何将这些数据绘成图形。（具体参看本文开头帖子中大伙儿的哀嚎。）正常情况下，要画这么多的点很花时间，还吃内存，更麻烦的是，还怎么画也不像Sam Derbyshire那幅“震撼”的图！要跨过这个槛，需要另一个技巧，好吧，其实这个技巧我是从这帖学来的：http://mathematica.stackexchange.com/questions/50839/smooth-peter-de-jong-attractor#comment149926_50839 （注意下面的答案最好也看看），这里的重点是，Sam那幅图中的颜色明暗到底表现的是什么？复平面上同一点处重合的点量？我们不难想像完全重合的点其实是很少的，是的，
颜色的明暗其实表现的只是某一个点附近小区域内的点量。如何决定这个小区域？对于这点，我们的计算机屏幕其实已经帮我们限制好了：
屏幕上的像素点是很有限的。于是这又反过来启发了我们：既然最终在屏幕上呈现的点数有限，我们又何必把所有的点都拿去画图？只要事先统计好每个像素内有多少点，再把这个统计结果拿去画图就好了嘛。哪个函数可以经济地画出这种原始数据点数极多的明暗图？我们有ArrayPlot。ArrayPlot可以调配色方案嗯。
另，最后画图时使用UnitStep剪掉了实轴上的点，因为影响视觉效果。
感觉其实还有些地方地方可以解说下不过现在也累了所以算了。
最后，其实即便抛开方程阶数，单论配色方案，我也没有完全“复现”人家的图。（小声：不过我自以为我选的方案更好看。）有谁有兴趣可以研究下怎么用Blend把原图的颜色混合出来。（应该不是某个Mathematica里现成的配色方案吧……）

赞！效果很好

很漂亮,用mathematica实现代码也很短
但是楼主只优化了时间复杂度, 没有考虑空间复杂度.
19次就要用2G内存的话,24次已然不现实了.

主体思想不变，局部小修改。这样就解决内存不够问题了，n=18时内存大约只需要几十M，运算时间基本不变。
Clear["Global`*"]
t0 = AbsoluteTime[];
modifiedroots[c_List] :=
Block[{a = DiagonalMatrix[ConstantArray[1., Length@c- 1], -1]},
a[[1]] = -c;
Eigenvalues[a]];
n = 18;
s = 1000;
data = ConstantArray[0., {s + 1, s + 1}];
For[i = 0, i < 2^n, i++,
data[[
#1, #
2]] += 1. & @@@
Round[1 +
s Function[z, {Im@z, Re@z}/4 + 0.5, Listable]@
modifiedroots[(-1)^IntegerDigits[i, 2, n]]];
]
AbsoluteTime[] - t0
ArrayPlot[UnitStep[80 - data] data, ColorFunction -> "AvocadoColors"]

稍微优化了空间复杂度，但是回复被系统删帖了，坐等人工恢复。

由于用For效率变低了点，所以加上了并行计算，下面的是双核。
Clear["Global`*"]
CloseKernels[]; LaunchKernels[2];
t0 = AbsoluteTime[];
ParallelTable[Module[{n, s, data, a, j},
n = 18; s = 1000; data = ConstantArray[0., {s + 1, s + 1}];
a = DiagonalMatrix[ConstantArray[1., n - 1], -1];
Do[
data[[
#1, #
2]] += 1. & @@@
Round[1 +
s Function[z, {Im@z, Re@z}/4 + 0.5, Listable]@
Eigenvalues[a[[1]] = (-1.)^IntegerDigits[j, 2, n]; a]
], {j, i - 1, 2^n - 1, 2}];
data
], {i, 2}];
l = Total@%;
AbsoluteTime[] - t0
ArrayPlot[UnitStep[80 - l] l, ColorFunction -> "AvocadoColors",
ImageSize -> 200]

那后来我也发现了这方法，不过为了得到可重复利用数据我选择的是导出数据。
也是2GB内存2GHz电脑做的20次(分4次导出) 。
导出导入相当耗时，并且当时用的是ListDensityPlot用了30min
代码太乱了，现在已经看不懂。。懒的折腾了。。
ListDensityPlot[image, Mesh -> False, Frame -> False,
ColorFunction -> (Blend[{{0, Black}, {1, Orange}, {60,
White}}, #] &)]

51分钟算出24次啦嘿嘿嘿~ 配置没变~
来个局部~

好强，小白飘过

图形怎么放大 看到细节啊

我看我还是贴代码吧, 并没有什么新东西, 就是结合了楼上各位的方法总结来的.
作为比较, 18次的不分块算了25秒, 分成2^10的块算了29秒. 内存方面之前说几个G是不分块的, 而且当然18, 19次并没有那么多, 用 MaxMemoryUsed 看的话, 18次280M左右, 19次也大概845M, 20次就多了, 1.75G左右. 所以上面说的就是分了18或19为一小块跑的24次.
roots[c_List] := Module[{a},
a = DiagonalMatrix[ConstantArray[1., Length[c] - 1], -1];
a[[All, -1]] = -c;
{-Im[
#], Re[#
]} & /@ Eigenvalues[a]] (* 负号是因为Image画图是上下颠倒的, 不过既然是对称就无所谓啦 *)
(* 不分块全部一起算 *)
data[pixels_, n_] := Module[{c, matrix, half, dp, max},
c = Tuples[{-1., 1.}, n];
half = 0.5 pixels; dp = pixels/3.6; (* Default range={{-1.8, 1.8},{-1.8,1.8}} *)
matrix = SparseArray[
Rule @@@ Tally[ (* 这里经 @苹果_nb 提示, 测试了下用 Rule @@ Transpose @ Tally , 次数低时要略快一点点, 次数高比如20次以上就不明显了 *)
Join @@ ParallelTable[1 + Floor[half + dp Select[roots[c[[i]]], #[[1]]^2 > 0 &]], {i, 2^n}, Method -> "CoarsestGrained"]
], {pixels, pixels}];
max = Max[matrix];
matrix = matrix/max
]
(* data[600,19] => time ≈ 54s (4 kernels, 2.4GHz, Mem < 850M) *)
(* 分成 2^nmin 的块来算 *)
data1[pixels_, n_, nmin_] :=
Module[{tmp, c, matrix, ind, imax, indmax, half, dp, max},
imax = 2^nmin; indmax = 2^Max[n - nmin, 0]; (* n>nmin *)
half = 0.5 pixels; dp = pixels/3.6; (* Default range={{-1.8,1.8},{-1.8,1.8}} *
matrix = SparseArray[{{1, 1} -> 0}, {pixels, pixels}];
Do[If[Mod[ind - 1, 2] == 0,
PrintTemporary[N[100 (ind - 1)/indmax, 4], "%"]]; (* 这是打印进度, 可以无视... *)
c = Table[(-1.)^IntegerDigits[i, 2, n], {i, (ind - 1) imax + 1, ind imax}];
matrix = matrix + SparseArray[
Rule @@@ Tally[
Join @@ ParallelTable[1 + Floor[half + dp Select[roots[c[[i]]], #[[1]]^2 > 0 &]], {i, imax}, Method -> "CoarsestGrained"]
], {pixels, pixels}];
, {ind, indmax}];
max = Max[matrix];
matrix = matrix/max
]
(* data1[600,19,18] => time ≈ 56s (4 kernels, 2.4GHz, Mem < 300M) *)
(* 画图, 当然这个图什么也看不见, 我是直接存成黑白png, 然后再慢慢Colorize的 *)
Image[data[pixels, n], ImageSize -> All]
Image[data1[pixels, n, nmin], ImageSize -> All]
(* Image的好处很多, 拼图的时候就会发现Image没有其他Plot那种边框(比如只画一半再用ImagePad拼图的时候) *)

学习了

没看懂啊，怎么做啊，教教我啊

ls可以先看吧主说的那个画 de jong attractor 的贴, 还是很有用的, 我也是先从那个学的, 以前只在画分形三维图的时候想过手动画像素点这种方法, 真的很省资源