Generalized series approximation
In a previous post I presented a neat form for the series approximation coefficient updates for the quadratic Mandelbrot set. Note that the series for the derivative is redundant, as you can just take the derivative of the main series term by term. Recently superheal on fractalforums.org asked about series approximation for other powers. I explored a bit with SageMath (which is based on Python) and came up with this code (empirically derived formula highlighted):
def a(k):
return var("a_" + str(k))
m = 10
for p in range(2,5):
f(c, z) = z^p + c
print(f)
g(C, Z, c, z) = (f(C + c, Z + z) - f(C, Z)).expand().simplify()
print(g)
h0(C, Z, c) = sum( a(k) * c^k for k in range(1, m) ).series(c, m)
print(h0)
h1(C, Z, c) = g(C, Z, c, h0(C, Z, c)).expand().series(c, m)
print(h1)
def x(C, Z, c, k):
return h1(C, Z, c).coefficient(c^k)
def y(C, Z, c, k):
return ((1 if k == 1 else 0) +
sum([binomial(p, t) * Z^(p-t) *
sum([Permutations(qs).cardinality() *
product([a(q) for q in qs])
for qs in Partitions(k, length=t).list()])
for t in range(1, k + 1)]))
for k in range(1, m):
print((x(C, Z, c, k) - y(C, Z, c, k)).full_simplify(), " : ", a(k)," := ", x(C, Z, c, k))
print()Example output:
(c, z) |--> z^2 + c (C, Z, c, z) |--> 2*Z*z + z^2 + c (C, Z, c) |--> (a_1)*c + (a_2)*c^2 + (a_3)*c^3 + (a_4)*c^4 + (a_5)*c^5 + (a_6)*c^6 + (a_7)*c^7 + (a_8)*c^8 + (a_9)*c^9 + Order(c^10) (C, Z, c) |--> (2*Z*a_1 + 1)*c + (a_1^2 + 2*Z*a_2)*c^2 + (2*a_1*a_2 + 2*Z*a_3)*c^3 + (a_2^2 + 2*a_1*a_3 + 2*Z*a_4)*c^4 + (2*a_2*a_3 + 2*a_1*a_4 + 2*Z*a_5)*c^5 + (a_3^2 + 2*a_2*a_4 + 2*a_1*a_5 + 2*Z*a_6)*c^6 + (2*a_3*a_4 + 2*a_2*a_5 + 2*a_1*a_6 + 2*Z*a_7)*c^7 + (a_4^2 + 2*a_3*a_5 + 2*a_2*a_6 + 2*a_1*a_7 + 2*Z*a_8)*c^8 + (2*a_4*a_5 + 2*a_3*a_6 + 2*a_2*a_7 + 2*a_1*a_8 + 2*Z*a_9)*c^9 + Order(c^10) 0 : a_1 := 2*Z*a_1 + 1 0 : a_2 := a_1^2 + 2*Z*a_2 0 : a_3 := 2*a_1*a_2 + 2*Z*a_3 0 : a_4 := a_2^2 + 2*a_1*a_3 + 2*Z*a_4 0 : a_5 := 2*a_2*a_3 + 2*a_1*a_4 + 2*Z*a_5 0 : a_6 := a_3^2 + 2*a_2*a_4 + 2*a_1*a_5 + 2*Z*a_6 0 : a_7 := 2*a_3*a_4 + 2*a_2*a_5 + 2*a_1*a_6 + 2*Z*a_7 0 : a_8 := a_4^2 + 2*a_3*a_5 + 2*a_2*a_6 + 2*a_1*a_7 + 2*Z*a_8 0 : a_9 := 2*a_4*a_5 + 2*a_3*a_6 + 2*a_2*a_7 + 2*a_1*a_8 + 2*Z*a_9 (c, z) |--> z^3 + c (C, Z, c, z) |--> 3*Z^2*z + 3*Z*z^2 + z^3 + c (C, Z, c) |--> (a_1)*c + (a_2)*c^2 + (a_3)*c^3 + (a_4)*c^4 + (a_5)*c^5 + (a_6)*c^6 + (a_7)*c^7 + (a_8)*c^8 + (a_9)*c^9 + Order(c^10) (C, Z, c) |--> (3*Z^2*a_1 + 1)*c + (3*Z*a_1^2 + 3*Z^2*a_2)*c^2 + (a_1^3 + 6*Z*a_1*a_2 + 3*Z^2*a_3)*c^3 + (3*a_1^2*a_2 + 3*Z*a_2^2 + 6*Z*a_1*a_3 + 3*Z^2*a_4)*c^4 + (3*a_1*a_2^2 + 3*a_1^2*a_3 + 6*Z*a_2*a_3 + 6*Z*a_1*a_4 + 3*Z^2*a_5)*c^5 + (a_2^3 + 6*a_1*a_2*a_3 + 3*Z*a_3^2 + 3*a_1^2*a_4 + 6*Z*a_2*a_4 + 6*Z*a_1*a_5 + 3*Z^2*a_6)*c^6 + (3*a_2^2*a_3 + 3*a_1*a_3^2 + 6*a_1*a_2*a_4 + 6*Z*a_3*a_4 + 3*a_1^2*a_5 + 6*Z*a_2*a_5 + 6*Z*a_1*a_6 + 3*Z^2*a_7)*c^7 + (3*a_2*a_3^2 + 3*a_2^2*a_4 + 6*a_1*a_3*a_4 + 3*Z*a_4^2 + 6*a_1*a_2*a_5 + 6*Z*a_3*a_5 + 3*a_1^2*a_6 + 6*Z*a_2*a_6 + 6*Z*a_1*a_7 + 3*Z^2*a_8)*c^8 + (a_3^3 + 6*a_2*a_3*a_4 + 3*a_1*a_4^2 + 3*a_2^2*a_5 + 6*a_1*a_3*a_5 + 6*Z*a_4*a_5 + 6*a_1*a_2*a_6 + 6*Z*a_3*a_6 + 3*a_1^2*a_7 + 6*Z*a_2*a_7 + 6*Z*a_1*a_8 + 3*Z^2*a_9)*c^9 + Order(c^10) 0 : a_1 := 3*Z^2*a_1 + 1 0 : a_2 := 3*Z*a_1^2 + 3*Z^2*a_2 0 : a_3 := a_1^3 + 6*Z*a_1*a_2 + 3*Z^2*a_3 0 : a_4 := 3*a_1^2*a_2 + 3*Z*a_2^2 + 6*Z*a_1*a_3 + 3*Z^2*a_4 0 : a_5 := 3*a_1*a_2^2 + 3*a_1^2*a_3 + 6*Z*a_2*a_3 + 6*Z*a_1*a_4 + 3*Z^2*a_5 0 : a_6 := a_2^3 + 6*a_1*a_2*a_3 + 3*Z*a_3^2 + 3*a_1^2*a_4 + 6*Z*a_2*a_4 + 6*Z*a_1*a_5 + 3*Z^2*a_6 0 : a_7 := 3*a_2^2*a_3 + 3*a_1*a_3^2 + 6*a_1*a_2*a_4 + 6*Z*a_3*a_4 + 3*a_1^2*a_5 + 6*Z*a_2*a_5 + 6*Z*a_1*a_6 + 3*Z^2*a_7 0 : a_8 := 3*a_2*a_3^2 + 3*a_2^2*a_4 + 6*a_1*a_3*a_4 + 3*Z*a_4^2 + 6*a_1*a_2*a_5 + 6*Z*a_3*a_5 + 3*a_1^2*a_6 + 6*Z*a_2*a_6 + 6*Z*a_1*a_7 + 3*Z^2*a_8 0 : a_9 := a_3^3 + 6*a_2*a_3*a_4 + 3*a_1*a_4^2 + 3*a_2^2*a_5 + 6*a_1*a_3*a_5 + 6*Z*a_4*a_5 + 6*a_1*a_2*a_6 + 6*Z*a_3*a_6 + 3*a_1^2*a_7 + 6*Z*a_2*a_7 + 6*Z*a_1*a_8 + 3*Z^2*a_9 (c, z) |--> z^4 + c (C, Z, c, z) |--> 4*Z^3*z + 6*Z^2*z^2 + 4*Z*z^3 + z^4 + c (C, Z, c) |--> (a_1)*c + (a_2)*c^2 + (a_3)*c^3 + (a_4)*c^4 + (a_5)*c^5 + (a_6)*c^6 + (a_7)*c^7 + (a_8)*c^8 + (a_9)*c^9 + Order(c^10) (C, Z, c) |--> (4*Z^3*a_1 + 1)*c + (6*Z^2*a_1^2 + 4*Z^3*a_2)*c^2 + (4*Z*a_1^3 + 12*Z^2*a_1*a_2 + 4*Z^3*a_3)*c^3 + (a_1^4 + 12*Z*a_1^2*a_2 + 6*Z^2*a_2^2 + 12*Z^2*a_1*a_3 + 4*Z^3*a_4)*c^4 + (4*a_1^3*a_2 + 12*Z*a_1*a_2^2 + 12*Z*a_1^2*a_3 + 12*Z^2*a_2*a_3 + 12*Z^2*a_1*a_4 + 4*Z^3*a_5)*c^5 + (6*a_1^2*a_2^2 + 4*Z*a_2^3 + 4*a_1^3*a_3 + 24*Z*a_1*a_2*a_3 + 6*Z^2*a_3^2 + 12*Z*a_1^2*a_4 + 12*Z^2*a_2*a_4 + 12*Z^2*a_1*a_5 + 4*Z^3*a_6)*c^6 + (4*a_1*a_2^3 + 12*a_1^2*a_2*a_3 + 12*Z*a_2^2*a_3 + 12*Z*a_1*a_3^2 + 4*a_1^3*a_4 + 24*Z*a_1*a_2*a_4 + 12*Z^2*a_3*a_4 + 12*Z*a_1^2*a_5 + 12*Z^2*a_2*a_5 + 12*Z^2*a_1*a_6 + 4*Z^3*a_7)*c^7 + (a_2^4 + 12*a_1*a_2^2*a_3 + 6*a_1^2*a_3^2 + 12*Z*a_2*a_3^2 + 12*a_1^2*a_2*a_4 + 12*Z*a_2^2*a_4 + 24*Z*a_1*a_3*a_4 + 6*Z^2*a_4^2 + 4*a_1^3*a_5 + 24*Z*a_1*a_2*a_5 + 12*Z^2*a_3*a_5 + 12*Z*a_1^2*a_6 + 12*Z^2*a_2*a_6 + 12*Z^2*a_1*a_7 + 4*Z^3*a_8)*c^8 + (4*a_2^3*a_3 + 12*a_1*a_2*a_3^2 + 4*Z*a_3^3 + 12*a_1*a_2^2*a_4 + 12*a_1^2*a_3*a_4 + 24*Z*a_2*a_3*a_4 + 12*Z*a_1*a_4^2 + 12*a_1^2*a_2*a_5 + 12*Z*a_2^2*a_5 + 24*Z*a_1*a_3*a_5 + 12*Z^2*a_4*a_5 + 4*a_1^3*a_6 + 24*Z*a_1*a_2*a_6 + 12*Z^2*a_3*a_6 + 12*Z*a_1^2*a_7 + 12*Z^2*a_2*a_7 + 12*Z^2*a_1*a_8 + 4*Z^3*a_9)*c^9 + Order(c^10) 0 : a_1 := 4*Z^3*a_1 + 1 0 : a_2 := 6*Z^2*a_1^2 + 4*Z^3*a_2 0 : a_3 := 4*Z*a_1^3 + 12*Z^2*a_1*a_2 + 4*Z^3*a_3 0 : a_4 := a_1^4 + 12*Z*a_1^2*a_2 + 6*Z^2*a_2^2 + 12*Z^2*a_1*a_3 + 4*Z^3*a_4 0 : a_5 := 4*a_1^3*a_2 + 12*Z*a_1*a_2^2 + 12*Z*a_1^2*a_3 + 12*Z^2*a_2*a_3 + 12*Z^2*a_1*a_4 + 4*Z^3*a_5 0 : a_6 := 6*a_1^2*a_2^2 + 4*Z*a_2^3 + 4*a_1^3*a_3 + 24*Z*a_1*a_2*a_3 + 6*Z^2*a_3^2 + 12*Z*a_1^2*a_4 + 12*Z^2*a_2*a_4 + 12*Z^2*a_1*a_5 + 4*Z^3*a_6 0 : a_7 := 4*a_1*a_2^3 + 12*a_1^2*a_2*a_3 + 12*Z*a_2^2*a_3 + 12*Z*a_1*a_3^2 + 4*a_1^3*a_4 + 24*Z*a_1*a_2*a_4 + 12*Z^2*a_3*a_4 + 12*Z*a_1^2*a_5 + 12*Z^2*a_2*a_5 + 12*Z^2*a_1*a_6 + 4*Z^3*a_7 0 : a_8 := a_2^4 + 12*a_1*a_2^2*a_3 + 6*a_1^2*a_3^2 + 12*Z*a_2*a_3^2 + 12*a_1^2*a_2*a_4 + 12*Z*a_2^2*a_4 + 24*Z*a_1*a_3*a_4 + 6*Z^2*a_4^2 + 4*a_1^3*a_5 + 24*Z*a_1*a_2*a_5 + 12*Z^2*a_3*a_5 + 12*Z*a_1^2*a_6 + 12*Z^2*a_2*a_6 + 12*Z^2*a_1*a_7 + 4*Z^3*a_8 0 : a_9 := 4*a_2^3*a_3 + 12*a_1*a_2*a_3^2 + 4*Z*a_3^3 + 12*a_1*a_2^2*a_4 + 12*a_1^2*a_3*a_4 + 24*Z*a_2*a_3*a_4 + 12*Z*a_1*a_4^2 + 12*a_1^2*a_2*a_5 + 12*Z*a_2^2*a_5 + 24*Z*a_1*a_3*a_5 + 12*Z^2*a_4*a_5 + 4*a_1^3*a_6 + 24*Z*a_1*a_2*a_6 + 12*Z^2*a_3*a_6 + 12*Z*a_1^2*a_7 + 12*Z^2*a_2*a_7 + 12*Z^2*a_1*a_8 + 4*Z^3*a_9
You can try it online. The first column should be all 0 if the formula is correct, which it seems to be for all the cases I've tried. But it remains to be proven rigourously that it is correct for all terms of all powers.
Efficiently implementing the general series coefficient update formula would probably need a multi-stage process: first (one-time cost given power and number of terms) calculate tables of constants (binomials, partitions as (index, multiplicity) pairs, partition permutation cardinalities). Second stage (once per iteration) calculate tables of powers of series coefficient variables (only going as far as the highest needed multiplicity for each index), and powers of Z. Third stage (once per series coefficient variable) combine all the powers and constants. Final stage, add 1 to the first variable.
It should be possible to generate OpenCL code for this at runtime. The expressions for each variable are of very different sizes but bundling a_k with a_(m-k) might give a more uniform load per execution unit.