Help about implementing the perturbation theory

[CFJH]

Hello,

I'm reading very interested the topic about sft with the significant speedups by using the perturbation theory.

http://www.fractalforums.com/announcements-and-news/superfractalthing-arbitrary-precision-mandelbrot-set-rendering-in-java/

Now I want to implement this in my own plain fractal generator written in C. But I havn't really understand how the perturbation theory is working.
Reviewing the java-code of sft is difficult to me because I'm not a java developer and some used structures are unknown to me.

Can someone explain me, how is it working or give a abstract (programming language indepent) code fragment ?

Jürgen

[Kalles Fraktaler]

I recommend that you take the time and read these two threads through and the linked pdf:
http://www.fractalforums.com/announcements-and-news/superfractalthing-arbitrary-precision-mandelbrot-set-rendering-in-java/
http://www.fractalforums.com/announcements-and-news/kalles-fraktaler-2/

Even if the Java code might be difficult to digest for a C/C++ programmer all info needed is there.

[Gluecker]

the most significant speed boost will come from equation 1 (see sft_maths.pdf from the SuperFractalThing thread). Since the java code provided has additional tweaks, its not really easy to fully understand.
I will try to explain what I got so far:

- calculate one point (the reference X) near or inside your field of view at full precision and store the real and imaginary part of every calculation at lower precision. This point should be inside the M-set, or at least very close to it.

- compute the initial difference between the point Y you want to render and the reference (using a lower precision)

- iterate over eq1 using and the stored values for X. Its important to start with the correct index for X. (Eg. I stored the starting coordinates of X as X[0], my first calculation is then X[1], so my iteration counter is initially set to 1)

- as usual the iteration ends when the iteration counter hits the iteration limit. You also want to check if the distance gets bigger than a certain limit, since the assumption that "all the numbers are 'small'" was made. works fine for me, but I am not done testing yet.

- when gets too big, we can calculate back to the current value of and go on with the standart mandelbrot calculation. This also seems to work with lower precision, but again, not sure about that.

hope that helps.
cheers

[CFJH]

Hallo
in the past days I had a litte time to implememt something about the pertubation theory.
Now I've an (nearly) working function which is currently using the equation (1) only.

But my problem is the bailout-value, where the interartion exits. The standard für a mandelbrot is 4 but I've to use a value 64 to get an unbroken image in all cases.

the attached Images shows
test1s : calculated with the classical method
test2s : image with bailout of 64
test3s : broken image with standard bailout of 4

Code:

int tMandel::CalcPointPertub(int ix, int iy) {
double D0r,D0i; // Delta 0
double DNr,DNi; // Delta n
double Xr,Xi;   // aktuelles Xn aus dem Array
double Yr,Yi;   // berechneter Yn-Wert
double Ti,Tr;   // temporär
unsigned long pIter;
int Inside=1;

  mpf_set_d(a,ix-0.5*SizeX);
  mpf_mul(a,dx,a);
  D0r=mpf_get_d(a);

  mpf_set_d(a,iy-0.5*SizeY);
  mpf_mul(a,dy,a);
  D0i=mpf_get_d(a);
  // Reference point in the center
  // ∆n+1 = 2Xn ∆n + ∆n^2 + ∆0
  DNr=D0r;
  DNi=D0i;
  for (pIter=0;pIter<MaxIter;pIter++) {
    Xr=arx[pIter]; // reference values real part
    Xi=ary[pIter]; // reference values imaginary part
    
    // Calculation 2Xn ∆n
    Tr = 2.0*(Xr*DNr - Xi*DNi);
    Ti = 2.0*(Xr*DNi + Xi*DNr);
    // + ∆n^2
    Tr = Tr + (DNr*DNr - DNi*DNi);
    Ti = Ti + 2.0*DNr*DNi;
    // + ∆0
    Tr += D0r;
    Ti += D0i;
    
    DNr = Tr;
    DNi = Ti;
    
    // ∆n = Yn − Xn  -> Yn = ∆n - Xn
    Yr = DNr - Xr;
    Yi = DNi - Xi;
    
    //Abort
    if ((Yr*Yr+Yi*Yi)>64.0) {  // <----- problem !!
      Inside=0;
      break;
    }
    if (Aborted) {
      Calculation=FALSE;
      return -1;
    }
  }
  if (Inside)
    return 0;
  else
    return pIter;
}

Why is the image broken with standard bailout value ?
I can't find where is the problem in my code. Can somebody give me a hint ?

[knighty]

// ∆n = Yn − Xn  -> Yn = Xn + ∆n

[Gluecker]

yea, Y = X+delta.
other than that, higher bailouts dont hurt since the assumption "|Z|> 2 escapes" also implies |Z|> 8 does.

[knighty]

Any zooms? 

[CFJH]

Thanks for the hint,

but this solves the problem not completly. Now I can decrease the bailout value to 16 (instrad of 64) but the rendered image still looks different against an image rendered with standard method. A lower value as 16 produces a clipped image.
Is this normal or have I another bug ?

Jürgen

[knighty]

Don't know. If it's approximately the same as the 3rd picture then I would say there still is an error in Yn calculation.

[Kalles Fraktaler]

Thanks for the hint,

but this solves the problem not completly. Now I can decrease the bailout value to 16 (instrad of 64) but the rendered image still looks different against an image rendered with standard method. A lower value as 16 produces a clipped image.
Is this normal or have I another bug ?

Jürgen

No it is NOT normal and yes you have a bug...
Sorry but I cannot find any error in you code either.
Here is my bailout validation, which looks very much like yours:

Code:

Di=D0i;
Dr=D0r;
for(antal=0;antal<m_nMaxIter;antal++){
	Dnr = (m_db_dxr[antal]*Dr - m_db_dxi[antal]*Di)*2 + Dr*Dr - Di*Di + D0r;
	Dni = (m_db_dxr[antal]*Di + m_db_dxi[antal]*Dr + Dr*Di)*2 + D0i;
	yr=m_db_dxr[antal] + Dr;
	yi=m_db_dxi[antal] + Di;

	Di = Dni;
	Dr = Dnr;

	if(yr*yr+yi*yi > 4)
		break;
}

Are you trying to do some coloring with the value of Yn? Because I have noticed that the value of Yn is unstable. It is accurate for calculating the bailout iteration, but it is not usable to do further coloring with the value since it get too far away for some locations.

[CFJH]

Hello,
now I could find my problem in the calculation 

My mistake was, that I compute the Yn values one step too late

correct:

Code:

      Tr = 2.0*Xr*DNr - 2.0*Xi*DNi +     DNr*DNr - DNi*DNi + D0r;
      Ti = 2.0*Xr*DNi + 2.0*Xi*DNr + 2.0*DNr*DNi           + D0i;    
      Yr=DNr+Xr;
      Yi=DNi+Xi;
      DNr = Tr;
      DNi = Ti;

wrong:

Code:

      Tr = 2.0*Xr*DNr - 2.0*Xi*DNi +     DNr*DNr - DNi*DNi + D0r;
      Ti = 2.0*Xr*DNi + 2.0*Xi*DNr + 2.0*DNr*DNi           + D0i;    
      DNr = Tr;
      DNi = Ti;
      Yr=DNr+Xr;
      Yi=DNi+Xi;

Thanks to all for the help!

Jürgen

[laser blaster]

Wow, this perturbation theory looks brilliant! I'm still struggling to wrap my head around the math, but if I understand it correctly, it should enable us to only have to iterate one point, store it's values at each iteration, and then discretely solve for the values of all other points in that area, without having to iterate them. That's huge! Even better, it seems like it could be applied to gpu rendering. I might be misinterpreting it, though.

[cKleinhuis]

You still need to iterate, but single precision is sufficient then

[Kalles Fraktaler]

You still need to iterate, but single precision is sufficient then

No single precision is not sufficient, at least I was not able to succeed with that

[Gluecker]

At a certain magnification delta gets too small (if the reference is near the screen center) to be calculated with single precision (minValue for float: 1.18 × 10^−38).
One could try to find good references further away, but thats not trivial and I guess there are also other limits...