END OF AN ERA, FRACTALFORUMS.COM IS CONTINUED ON FRACTALFORUMS.ORG

I restored the image in my previous post (my internet connection died while I was uploading the the post so the image didn't get loaded.)

I guess with Mathematica you can find inverses to more of them than I previously thought.

From his description, I still think he's using a forward iteration method. It's a variation of Cauchy convergence, comparing iteration of z and z+epsilon. The trick is to map C to the Riemann sphere for the distance metric, so that points both inside and outside the set converge, while points sufficiently close to the set don't. Here's someone doing a similar thing but with more color:
http://de.wikipedia.org/wiki/Benutzer:Georg-Johann/Mathematik

Did you get the Glynn set working?

hmm, perhaps with a generalisation to 3d orbits (via quaternions?) we can finally get a distance estimator for arbitrary functions... looks like this would be the space in which the magic would happen: http://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Quaternion_rotation_operations

if the cauchy series diverges, then one can use as a measure of divergence, which could be used to bound the distance to the surface from .

I'm fairly sure it's the root checking. I was also getting variations on such circular patterns before I had Peter Liepa's test for root validity:
theta = theta + TWOPI*(branch - floor((theta+limit)/TWOPI)) ;
if (theta > -limit && theta < limit){
theta = theta/parmC ;
} else {
bNoRoot = TRUE ;
}
With that test I was able to get it working in standard chaos game, before even trying MIIM. I didn't recognize anything in the Mathematica code to handle that. No root just triggers backtracking.

Doing the math by hand I can only handle inverses up to cubics, unless they resolve nicely. I've implemented (z^n+c)/(z^n-c) and will try one or two cubics.


I don't see any reason why not. Before I found that link, I did a quick test of the idea, without transforming to the Riemann sphere, for 2D quadratic Julia. It produced a good filled set, although it needs oversampling for antialiasing. The transformation is needed to separate both inside and outside from the set.

Thank you Garth. It works! I didn't think that complex roots could ever be invalid so I ignored that part of the code.

I may be the only Mathematica user here, but here's my corrected Mathematica code for the sake of posterity:

Pow[z_, n_, k_] := Module[{theta = Arg[z]}, theta = n(theta + 2Pi (k - Floor[(theta/Pi +Abs[1/n])/2])); If[Abs[theta] > Pi, Null, Abs[z]^n Exp[I theta]]];
power = 1.5; nroot = 3; zc = -0.2; zlist = {}; dzmax = 25.0; imax =1000; z = Table[-0.61, {imax}]; dz = Table[1, {imax}]; roots = Table[1, {imax}]; i = 2;
While[i > 1, z[[ i]] = Pow[z[[ i - 1]] - zc, 1/power, roots[[ i]] - 1]; If[z[[ i]] === Null, prune = True, dz[[ i]] = Abs[power]Abs[z[[ i]]]^(power - 1)dz[[ i - 1]]; zlist = Append[zlist, z[[ i]]]; prune = (i == imax || dz[[ i]] > dzmax)]; If[prune, While[i > 1 && roots[[ i]] == nroot, roots[[ i]] = 1; i--]; roots[[ i]]++, i++; roots[[ i]] = 1]];
ListPlot[{Re[ #], Im[ #]} & /@ zlist, PlotStyle -> PointSize[0.005], AspectRatio ->Automatic, Axes -> None]

ahh, that does make sense! and explains the problems i had with my rational-exponent inverse-iteration back in 2006... i hope to revisit inverse iteration for polynomials sometime soon, perhaps this weekend. btw, i think you have excellent taste in choosing the conjugating (z^n + a) / (z^n + b) iteration as a basic transformation (or holon in xenodream parlance) 

thanks again guys for all the great info in this thread that cauchy series approach is really cool and i hope to try it out before too long!

OK, well here are my first attempts to create 3D Glynn fractals using the triplex formulas. The first image is a slice that was ray traced using the escape time method. But it doesn't show any tree branch details so it's pretty boring. The other images are point clouds that were calculated using the MIIM method. Finding all the unique valid roots was a challenge. In general, I found that there are at least 2 roots when x < 0 and at least one root when x > 0. But there is also a strange conic section region containing 3 roots when z² > x² + y². It should also be noted that convergence is very sensitive to the initial seed value. I found the best results when the seed value is inside one of the cones of the conic section. But this only finds one half of the overall fractal. You can find the other half by placing the seed in the other cone, but I preferred to leave this part out, because it gives a nice cross-section so you can see inside the fractal. In these images c={-0.2,0,0}.

Here is some Mathematica code:
(* runtime: 15 seconds *)
power = 1.5; norm[x_] := x.x;
TriplexPow[{x_, y_, z_}, n_] := Module[{r = Sqrt[x^2 + y^2 + z^2], theta = n ArcTan[x, y], phi},phi = n ArcSin[z/r]; r^n{Cos[theta] Cos[phi], Sin[theta]Cos[phi], Sin[phi]}];
NumRootsGlynn[{x_, y_, z_}] := If[z^2 > x^2 + y^2, 3, If[x < 0, 2, 1]];
TriplexRootGlynn[{x_, y_, z_}, k_] := Module[{r = Sqrt[x^2 + y^2 + z^2], ktheta, kphi, theta, phi}, {ktheta, kphi} = If[k == 0, {0, 0}, If[z^2 > x^2 + y^2, If[k == 1, If[x <0 && y < 0, {2, 0}, {0.5, If[z > 0, 0.5, 2.5]}], If[x < 0 && y > 0, {1, 0}, {2.5, If[z > 0, 0.5, 2.5]}]], {If[y < 0,
2, 1], 0}]]; theta = (ArcTan[x, y] + ktheta Pi)/power; phi = (ArcSin[z/r] +kphi Pi)/power; r^(1/power){Cos[theta] Cos[phi], Sin[theta]Cos[phi], Sin[phi]}];
c = {-0.2, 0, 0}; imax = 100; dpmax = 50.0; p = Table[{-0.61, 0.0, -0.42}, {imax}]; dp = Table[1, {imax}]; roots = Table[0, {imax}]; plist = {}; i = 2;
While[i > 1 && Length[plist] < 10000, p[[ i]] = TriplexRootGlynn[p[[ i - 1]] - c, roots[[ i]]]; dp[[ i]] = Abs[power]Sqrt[norm[TriplexPow[p[[ i]], power - 1]]]dp[[ i - 1]];plist = Append[plist, p[[ i]]]; prune = (i == imax || dp[[ i]] > dpmax); If[prune, While[i > 1 && roots[[ i]] == NumRootsGlynn[p[[ i - 1]] - c] - 1, roots[[ i]] = 0; i--]; roots[[ i]]++, i++; roots[[ i]] = 0]];
Show[Graphics3D[{PointSize[0.005], Point /@ plist}]]

For those of you who don't use Mathematica, here is some C++ code:
//Glynn roots: {x,y,z}^(1/1.5) = {x,y,z}^(2/3)
int NumRootsGlynn(triplex p) { // finds number of branches for a given triplex p
 bool cone=(sqr(p.z)>sqr(p.x)+sqr(p.y));
 return (cone?3:(p.x<0.0?2:1));
}
triplex RootGlynn(triplex p, int k) {double n=1.5;
 bool cone=(sqr(p.z)>sqr(p.x)+sqr(p.y)); double ktheta, kphi;
 if(k==0) {ktheta=kphi=0.0;} else {
  if(!cone) {ktheta=(p.y<0.0?2.0:1.0); kphi=0.0;} else {
   if(k==1) {
    if(p.x<0.0 && p.y<0.0) {ktheta=2.0; kphi=0.0;} else {ktheta=0.5; kphi=(p.z>0.0?0.5:2.5);}
   }
   else {
    if(p.x<0.0 && p.y>0.0) {ktheta=1.0; kphi=0.0;} else {ktheta=2.5; kphi=(p.z>0.0?0.5:2.5);}
   }
  }
 }
 double r=abs(p), theta=(atan2(p.y,p.x)+ktheta*pi)/n, phi=(asin(p.z/r)+kphi*pi)/n;
 double cosphi=cos(phi);
 return pow(r,1.0/n)*triplex(cos(theta)*cosphi,sin(theta)*cosphi,sin(phi));
}