Welcome to Fractal Forums

Hello

Found a simple way to compute DE for trig. mandelbulbs.
I think it is new since I have not seen it on the forums?
I mention it anyway.

It is based on:

Dist = 0.5*|Zn|*log(|Zn|)/|DZn|

and

Zn+1 = Zn^P +C

DZn+1 = P*Zn^(P-1)*DZn +1

Thus nothing new so far!!

What is new (i think?) is the way to compute the product P*Zn^(P-1)*DZn

With Zn+1 defined as follows (wikipedia like def. of spherical coordinates, works with usual def. to.)
-----------------------------------------------------------------------------------------------
r=length(Zn)
t=acosf(Zn.z/r)
p=atan2f(Zn.y,Zn.x)

st = (r^P)*sin(P*t)

Zn+1.x = st*cos(P*p)+C.x
Zn+1.y = st*sin(P*p)+C.y
Zn+1.z = cos(P*t)*r^P+C.z

Now the derivative DZn+1 is:

l = length(DZn)

t = t *(P-1) + acos(DZn.z/l)
p = p*(P-1) + atan2(DZn.y,DZn.x)

rd = l*P*r^(P-1)

st = rd*sin(t)

DZn+1.x = st*cos(p)+1
DZn+1.y = st*sin(p)
DZn+1.z = cos(t)*rd


Thus I add the angles of Zn-1^P and DZn together and multiply their lengths..
to create the product vector. Exactly the  same way you can do for complex vectors.

I have tested it for the Z^2 mandel bulb and it works  ;D

I had to set:

DZ0 = (0,0,1)

To  get correct results.. if you use different def. of spherical coords you might need to try (1,0,0) and (0,1,0).

Also I set Z0 = C.

Can not show results yet, cause I am on a 44kbit modem  :( at the moment.

That's essentialy what's being used ;)

You should also note that because there is no "+1" involved in the derivative for Julia Sets of the form a*z^p+c then all you actually need is the running product of the magnitude i.e. you just do:

Dr = Dr*length(Zn)

where initial Dr is 1.0

Of course that optimisation won't adapt for exp(z)+c or z^p+z^p1+c etc.

Ok..so it was not new  :( Took me a week to figure it out..  :)
You say essentially..what is done different?

Well I'm pretty sure I'm using a start for the running derivative of (1,0,0) in all cases and assuming you're starting with initial z as (0,0,0) then that's correct since on the first iteration for the Mandelbrot you get 0 + 1 for the running derivative anyway.

I see.
Im not sure why I need different start values..but it looks alright in the end.

Why did you make that binary search algorithm for DE when you can use this method?
Is binary search better ?

I think you're confusing my delta DE method and using a binary search when solid is hit ?

I created the delta DE method for formulas other than the Mandelbulb before Paul Nylander extended Daniel White's original formula to higher powers finding the Mandelbulb - the other formulas concerned, although simple, are not so easy to get the analytical derivatives for - in fact for the Mandelbrot version they really need the full Jacobian which is both complicated and slow to calculate.
I used the delta DE method for my earlier renders of the Mandelbulb because it worked and because I'm a little lazy when it comes to maths (rather than programming) and I waited for someone else to come up with an analytical DE for the Mandelbulb using the derivative (I'm not sure who came up with that first but Jos Leys described it to me).

As for the binary search - that's a separate subject, when you use DE (delta or derivative) but you set the solid threshold to a high value, i.e. for a low resolution render, then without using a binary search (by initially stepping backwards 1/2 the last step when solid is hit) you get "stepping" in the render.

I should add that for some renders of the Mandelbulb using the latest optimised version of my Ultra Fractal formula (not released yet) then in some cases using delta DE is actually faster than using derivative DE - but not when using the magnitude only method for the derivative.

I should add that the delta DE method will work with *any* formula that you can get accurate smooth iteration values for.

Hi

Because analytic DE method is limited to fractals where we can find analytic derivatives, yesterday I tried to fight with delta DE methods. David, I tried to analyse your delta DE method but I can't understand some issues (for example I don't know how to calculate smooth iteration value). I decided to find method by myself.
In my method I calculate central difference instead of analytic derivative of z value, using z value for two point in very close distance (1e-10). I will describe my method in following pseudo-code:

delta = 1e-10
z1 = CalculateIterations(a,b,c)            //z1, z2 last values of z (after last iteration)
z2 = CalculateIterations(a+delta,b,c)   // I tried also (a+delta, b+delta, c+delta)
r = |(z1+z2)/2|
dr = |z2 - z1| / delta
distance = r * log(r) / dr

I tried this method in my 3D rendered. It works but there was lots of artifacts. I attached result ("delta DE test.jpg")

I made some investigations and tests on 2D slices. This method works very fine only for plane z = 0. When I tried to render another slices, for example z = 0.1 there was artifacts. I fought that I made some mistake in program but I tested it on quaternionic formula. For quaternions it works fine in full range of z.
I attached results of tests. On left side there are slices for z = 0 and on right side for z = 0.1.

David, you have very big experience with delta DE method. Do you know what I'm doing wrong?

I think so, you need to use 2 points *on the ray* i.e. at vp + alpha*d (the current step position) and vp + (alpha+1e-10)*d because the method you're using (like mine) is directional unlike the analytical method which gives the minimum distance in any direction.

(vp = viewpoint (ray start), alpha = step, d = ray's unit direction vector)

As to the smooth iteration value, use:

smooth = final iter + 1 + (log(log(bailout))-log(log(final magnitude)))/log(power)

Or (for much more accurate results) instead of log(power) use log(log(|z|)/log(|zold|)) when final iter is >1 (where zold is the penultimate value of z i.e. z from the iteration prior to bailout).

Actually for delta DE purposes because the difference of the 2 smooth iteration values is actually used then when using /log(power) then both the +1 term and the log(log(bailout)) term can be omitted - if using the log(log(|z|)/log(|zold|)) method then the +1 can be omitted but not the bailout term. I don't omit anything in my code because the overhead is small anyway and often the full correct iteration value will be needed for colouring.
Note that the +1 term is only required anyway if bailout on the first iteration results in a final iter of zero rather than 1.

I begin to understand the history..

Why can you just not step back a half step and set the threshold to a lower value.. ?
I guess I dont understand fully..that would be the same as having a low thresold to begin with or..

Binary search is for finding exact place where distance equals to DE threshold. Last DE step is always little too close (distance is less than DE threshold) and you have to cancel last step and go half step forward and check distance. When distance is to high you have to go forward or when distance is too small you have to go back again. In my programs I'm ending searching when distance is in range 90% - 100% of DE threshold. It is completely enough for all shading algorithms.

Thanks David for reply.
In my 3D renderer I used direction vector for delta. I also tested different directions of delta vector on slices but it always looks nearly the same. Very interesting for me is that it works properly for z = 0 and quaternions. I think our 3D/4D number systems don't have authentic derivatives and this a reason of artefacts.  Only complex numbers and quaternions has derivatives - analytical and accurate numerical.

I will try with smooth iteration values according to your description.

Yes, that would be the same as having a low threshold to begin with :)

Sometimes setting the DE large enough to avoid aliasing causes the stepping i.e. setting the DE threshold too small on a lower resolution image can cause aliasing but if you set it too high then you can get stepping.

Here's a gross example of the difference on a DE threshold animated quaternion (note that with my latest code the jagged artifacts on the version using the binary search have gone altogether):

http://makinmagic.deviantart.com/art/Plain-DE-animation-127361800 (http://makinmagic.deviantart.com/art/Plain-DE-animation-127361800)

http://makinmagic.deviantart.com/art/Using-DE-Binary-Search-127362040 (http://makinmagic.deviantart.com/art/Using-DE-Binary-Search-127362040)

Have you implimented the array to check that the new step distance is not greater than the current minimum for that iteration count on the ray ?
Without that you'll always get gross overstepping in some places when using the delta method.

No, I didn't implemented it. But I think for it I have to change final condition if I implement it. Now I'm searching fractal boundary until I achieve proper distance. Maybe will it be better when I finish searching when I find required value of iteration count? Now when my program meets some place where there is wrong value of distance (for example iteration count in this place is very low) and this is unfortunately lees than distance threshold, it finishes searching. When DE values are accurate everything is very simple but for power 2 fractals unfortunately are very inaccurate, also analytical.
Tomorrow I will try also with smooth iteration value instead of last value of z.

I see.
What about setting a scale factor for DE that underestimates the distance more?
Or will it still overestimate the last step if you scale ?

Using a high "accuracy" setting i.e. actually stepping distances much smaller than the actual distance estimate effectively performs the same function as the binary search, but to be as accurate at finding the exact threshold surface as a binary search the step distances would have to be so much smaller than the distance estimates as to make rendering incredibly slow.
For instance let's say it takes 20 steps to hit the surface using say distance estimate/2 step distances, if you then perform 4 binary steps to find the exact surface (24 steps in total) then your maximum distance from the threshold is distance estimate/2*(1/2)^4 i.e. distance estimate/32, so to get the same accuracy by reducing the step sizes you'd have to use step sizes of distance estimate/32 instead of distance estimate/2 so instead of 20 steps to hit the surface it could easily take over 200.

I made investigation why Delta DE causes artifacts on images. I found that artefacts was in areas where last value of theta (theta = atan2(y,x)) is changing value very fast. It caused very high values of dr (dr = |z2 - z1| / delta) because there was very big difference between z2 and z1 (nearly the same absolute value but almost opposite direction). I attached images with artefacts and theta value.

Regarding to this I changed a little algorithm for calculating Delta DE value:

delta = 1e-10
z1 = CalculateIterations(a,b,c)            //z1, z2 last values of z (after last iteration)
z2 = CalculateIterations(a+delta,b,c)   // I tried also (a+delta, b+delta, c+delta)
r = |(z1+z2)/2|
dr = (|z2| - |z1|) / delta
distance = r * log(r) / dr

Now dr is calculated as difference between absolute values. Results was better but estimator was very "directional". Without shading fractal looks well but when I tried to use shading algorithms, fractal has very unnatural surface (fractal was deformed along axis perpendicular to screen plane).

Finally I decided to calculate distance in three directions and calculate resultant distance dr = sqrt(dr.x^2 + dr.y^2 + dr.z^2). It needs to calculate iterations for 4 points but result is very accurate. Images are without any artifacts when programs don't have any additional procedures for avoiding errors. There are some examples bellow. For all images I used exactly the same Delta DE algorithm. IMHO I can say that this is very universal algorithm  :D.

Hypercomplex formula:
(http://www.fractalforums.com/gallery/1/640_06_01_10_8_37_52.jpg)

Trigonometric formula:
(http://www.fractalforums.com/gallery/1/640_06_01_10_8_38_34.jpg)

Quaternionic formula:
(http://www.fractalforums.com/gallery/1/640_06_01_10_8_39_22.jpg)

very nice renderings!

Excellent ! I think that's going to be faster than my smooth iteration delta DE - though it would still need testing on something more transcendental like quaternionic sin(z)+c for example or something including conditionals such as 3D+ versions of the Barnsley formulas :)
Have you tried it out on complex Julibrots ?
Also please check it out on disconnected Julias, or at least Julia Sets where back areas are visible through holes in front areas - you may find an issue with over-stepping (which using the min distance on each iteration array method will fix - you should note that the array method even fixes such problems that occur when using analytical DE and the overhead is considerably less than you might think - at least on a CPU, and is considerably more optimum than reducing all step distances to a point where such errors go away).

@Buddhi: Did you have any overstepping problems with the Julias you recently posted to the gallery ? (i.e. did you need to reduce the step distances more than normal)
I'm going to try your method over the weekend anyway - also watch out for a post from me (in "Mandelbulb Implimentation") describing a good position/framing of the degree 8 "-sine" Mandelbulb for really testing out rendering algorithms.

@David Makin: Unfortunately I had overstepping problems and I reduced step distances. The highest overstepping is on those places where fractal is very flattened. I'm interpreting this as the space around those places is also flattened and dimensions are deformed.
Plus of this method is that I didn't observed any artifacts on images and I don't have to use trigonometric versions of formulas, which are very slow.
I'm thinking now about some self regulating factor for reducing distance. If I have some results I will share.

Example of using Delta DE method for Julia z^2+c fractals. Factor for reducing step distance was 0.2. It was not so fast to render but any flat details was misses.

(http://www.fractalforums.com/gallery/1/640_08_01_10_10_34_42.jpg)

Buddhi, I'm trying to reproduce your image in UF.
Can you tell me what Julia seed you used exactly?

c = {0.51766599571757, -0.27872405784328, 0}

Thanks!

here is the same Julia using Analytical DE in UF. I must say yours looks a lot better!
How about timing? This 640*640 image took 35 seconds to render onscreen (so no anti-aliasing)

I plan to test your method for DE also.
Just to be clear : you iterate a point Px,Py,Pz which gives you an R after bailout. (R=sqrt(Px^2+Py^2+Pz^2) )
You then iterate also <Px+delta,Py,Pz>, <Px,Py+delta,Pz>, <Px,Py,Pz+delta>, which gives you Rx,Ry, Rz as values after bailout..
is then drx=|R-Rx|/delta, dry=|R-Ry|/delta, drz=|R-Rz|/delta, to give you dR=sqrt(drx*drx+dry*dry+drz*drz) 
and DE=R*log(R)/dR..?
At first sight this would not give a good DE..

Or do I have it totally wrong?

Using my program it takes 20 seconds to render 640x640 image (step distance was 0.3*DE).


As I write before I changed way to calculate differences, because numeric errors was to high on areas where angle of triplex number was changing very fast. It should be:

(EDITED! There was mistake in formulas)
drx=|R-Rx|/delta
dry=|R-Ry|/delta
drz=|R-Rz|/delta

Rest of formulas looks properly.

My image looks very good because I used very small factor for step distance (step = 0.1*DE) and It was rendered in 2560x2560 resolution with ambient occlusion based on many rays. Render time was 1 hour.

Hi guys, just to point out that if you're comparing timings you really need to state the computing power used ;)

OK, and thanks :)

I'm not sure I understand, Buddhi, R being a distance from the origin, is this not always positive?

If course  :o. My mistake. I don't know why I thought that R is the last value of triplex number after iterations. It should be:
drx=|R-Rx|/delta
dry=|R-Ry|/delta
drz=|R-Rz|/delta
It is exactly the same as you wrote. It means that your formulas are completely the same as mine.

double th = atan2(y, x);
double ph = -atan2(z, sqrt(x * x + y * y));
x = r2 * cos(p * ph) * cos(p * th) + par.a;
y = r2 * cos(p * ph) * sin(p * th) + par.b;
z = r2 * sin(p * ph) + par.c;

Managed to get a similar view of the same Julia...

(http://www.fractalforums.com/gallery/1/141_09_01_10_5_40_42.jpg)

Rendered within UF @640*640 in 27secs on my 3GHz P4HT using the current optimised version of my formula, just to put that in context exactly the same image rendered using full trig for both the iterate and the derivative takes nearly 5 minutes.
Rendered using analytical DE, threshold 1e-3, step sizes = 0.45*DE.

Here's part of the optimised code:

              magn = sqrt(magn)
              dr = @mpwr*(ph=sqrt((th=|dzri|) + sqr(dzj)))
              if (r=|zri|)>0.0
                zjk = ((r=sqrt(r)) + flip(zj))/magn
                dr = dr*magn^(@mpwr-1.0)
                zri = zri/r
                if @mpwr==round(@mpwr) && abs(@mpwr-1.0)<32
                  if @mpwr==1.0
                    zri1 = (1,0)
                    zjk1 = (1,0)
                  elseif @mpwr==0.0
                    zri1 = 1.0/zri
                    zjk1 = 1.0/zjk
                  else
                    if @mpwr-1.0<0.0
                      zri1 = 1.0/zri
                      zjk1 = 1.0/zjk
                    else
                      zri1 = zri
                      zjk1 = zjk
                    endif
                    if abs(@mpwr-1.0)==2.0
                      zri1 = sqr(zri1)
                      zjk1 = sqr(zjk1)
                    elseif abs(@mpwr-1.0)==3.0
                      zri1 = zri1^3
                      zjk1 = zjk1^3
                    elseif abs(@mpwr-1.0)==4.0
                      zri1 = zri1^4
                      zjk1 = zjk1^4
                    else
                      if (abs(@mpwr-1.0)%4)>1
                        zp2 = sqr(zri1)
                        z1p2 = sqr(zjk1)
                      endif
                      if (abs(@mpwr-1)%8)>3
                        if (abs(@mpwr-1)%4)>1
                          zp4 = sqr(zp2)
                          z1p4 = sqr(z1p2)
                        else
                          zp4 = zri1^4
                          z1p4 = zjk1^4
                        endif
                      endif
                      if (abs(@mpwr-1)%16)>7
                        if (abs(@mpwr-1)%8)>3
                          zp8 = sqr(zp4)
                          z1p8 = sqr(z1p4)
                        elseif (abs(@mpwr-1)%4)>1
                          zp8 = zp2^4
                          z1p8 = z1p2^4
                        else
                          zp8 = sqr(zri1^4)
                          z1p8 = sqr(zjk1^4)
                        endif
                      endif
                      if abs(@mpwr-1)>15
                        if (abs(@mpwr-1)%16)>7
                          zp16 = sqr(zp8)
                          z1p16 = sqr(z1p8)
                        elseif (abs(@mpwr-1)%8)>3
                          zp16 = zp4^4
                          z1p16 = z1p4^4
                        elseif (abs(@mpwr-1)%4)>1
                          zp16 = sqr(zp2^4)
                          z1p16 = sqr(z1p2^4)
                        else
                          zp16 = (zri1^4)^4
                          z1p16 = (zjk1^4)^4
                        endif
                      endif
                      if abs(@mpwr-1)==31
                        zri1 = zp16*zp8*zp4*zp2*zri1
                        zjk1 = z1p16*z1p8*z1p4*z1p2*zjk1
                      elseif abs(@mpwr-1)==30
                        zri1 = zp16*zp8*zp4*zp2
                        zjk1 = z1p16*z1p8*z1p4*z1p2
                      elseif abs(@mpwr-1)==29
                        zri1 = zp16*zp8*zp4*zri1
                        zjk1 = z1p16*z1p8*z1p4*zjk1
                      elseif abs(@mpwr-1)==28
                        zri1 = zp16*zp8*zp4
                        zjk1 = z1p16*z1p8*z1p4
                      elseif abs(@mpwr-1)==27
                        zri1 = zp16*zp8*zp2*zri1
                        zjk1 = z1p16*z1p8*z1p2*zjk1
                      elseif abs(@mpwr-1)==26
                        zri1 = zp16*zp8*zp2
                        zjk1 = z1p16*z1p8*z1p2
                      elseif abs(@mpwr-1)==25
                        zri1 = zp16*zp8*zri1
                        zjk1 = z1p16*z1p8*zjk1
                      elseif abs(@mpwr-1)==24
                        zri1 = zp16*zp8
                        zjk1 = z1p16*z1p8
                      elseif abs(@mpwr-1)==23
                        zri1 = zp16*zp4*zp2*zri1
                        zjk1 = z1p16*z1p4*z1p2*zjk1
                      elseif abs(@mpwr-1)==22
                        zri1 = zp16*zp4*zp2
                        zjk1 = z1p16*z1p4*z1p2
                      elseif abs(@mpwr-1)==21
                        zri1 = zp16*zp4*zri1
                        zjk1 = z1p16*z1p4*zjk1
                      elseif abs(@mpwr-1)==20
                        zri1 = zp16*zp4
                        zjk1 = z1p16*z1p4
                      elseif abs(@mpwr-1)==19
                        zri1 = zp16*zp2*zri1
                        zjk1 = z1p16*z1p2*zjk1
                      elseif abs(@mpwr-1)==18
                        zri1 = zp16*zp2
                        zjk1 = z1p16*z1p2
                      elseif abs(@mpwr-1)==17
                        zri1 = zp16*zri1
                        zjk1 = z1p16*zjk1
                      elseif abs(@mpwr-1)==16
                        zri1 = zp16
                        zjk1 = z1p16
                      elseif abs(@mpwr-1)==15
                        zri1 = zp8*zp4*zp2*zri1
                        zjk1 = z1p8*z1p4*z1p2*zjk1
                      elseif abs(@mpwr-1)==14
                        zri1 = zp8*zp4*zp2
                        zjk1 = z1p8*z1p4*z1p2
                      elseif abs(@mpwr-1)==13
                        zri1 = zp8*zp4*zri1
                        zjk1 = z1p8*z1p4*zjk1
                      elseif abs(@mpwr-1)==12
                        zri1 = zp8*zp4
                        zjk1 = z1p8*z1p4
                      elseif abs(@mpwr-1)==11
                        zri1 = zp8*zp2*zri1
                        zjk1 = z1p8*z1p2*zjk1
                      elseif abs(@mpwr-1)==10
                        zri1 = zp8*zp2
                        zjk1 = z1p8*z1p2
                      elseif abs(@mpwr-1)==9
                        zri1 = zp8*zri1
                        zjk1 = z1p8*zjk1
                      elseif abs(@mpwr-1)==8
                        zri1 = zp8
                        zjk1 = z1p8
                      elseif abs(@mpwr-1)==7
                        zri1 = zp4*zp2*zri1
                        zjk1 = z1p4*z1p2*zjk1
                      elseif abs(@mpwr-1)==6
                        zri1 = zp4*zp2
                        zjk1 = z1p4*z1p2
                      elseif abs(@mpwr-1)==5
                        zri1 = zp4*zri1
                        zjk1 = z1p4*zjk1
                      elseif abs(@mpwr-1)==4
                        zri1 = zp4
                        zjk1 = z1p4
                      elseif abs(@mpwr-1)==3
                        zri1 = zp2*zri1
                        zjk1 = z1p2*zjk1
                      elseif abs(@mpwr-1)==2
                        zri1 = zp2
                        zjk1 = z1p2
                      endif
                    endif
                  endif
                else
                  zri1 = zri^(@mpwr-1.0)
                  zjk1 = zjk^(@mpwr-1.0)
                endif
                r = sqrt(th)
                zri = zri*zri1
                dzjk = r + flip(dzj)
                zjk = zjk*zjk1
                if r>0.0
                  dzri = zri1*dzri/r
                else
                  dzri = zri1
                endif
                if ph>0.0
                  dzjk = zjk1*dzjk/ph
                else
                  dzjk = zjk1
                endif
                magn = magn^@mpwr
                dzri = dr*real(dzjk)*dzri
                if @plus
                  dzj = dr*imag(dzjk)
                else
                  dzj = -dr*imag(dzjk)
                endif
                if @fractaltype==1
                  dzri = dzri + 1.0
                endif
                if @plus
                  zj = cj + magn*imag(zjk)
                else
                  zj = cj - magn*imag(zjk)
                endif
                zri = magn*real(zjk)*zri + cri

Thanks - I realised it was just a matter of changing viewing angle and perspective to get the image so I removed the post asking which formula it was ;)

Apologies - that code's only used for the Mandelbrot !
Here's part of the Julia version:

            elseif (r=|zri|)>0.0
              magn = sqrt(magn)
              if @usede
                dr = @mpwr*dr*magn^(@mpwr-1.0)
                if @fractaltype==1
                  dr = dr + 1.0
                endif
              endif
              if @mpwr==0.0
                zri = 1.0 + cri
                zj = cj
              elseif @mpwr==1.0
                zri = zri + cri
                if @plus
                  zj = cj + zj
                else
                  zj = cj - zj
                endif
              elseif @mpwr==2.0
                zri = sqr(zri)*(1.0 - sqr(zj)/r) + cri
                if @plus
                  zj = cj + 2.0*zj*sqrt(r)
                else
                  zj = cj - 2.0*zj*sqrt(r)
                endif
              elseif @mpwr==3.0
                zri = (zri^3)*(1.0 - 3.0*sqr(zj)/r) + cri
                if @plus
                  zj = cj - zj^3 + 3.0*zj*r
                else
                  zj = zj^3 - 3.0*zj*r + cj
                endif
              elseif @mpwr==4.0
                ztemp = ((r=sqrt(r)) + flip(zj))^4
                zri = real(ztemp)*(zri/r)^4 + cri
                if @plus
                  zj = cj + imag(ztemp)
                else
                  zj = cj - imag(ztemp)
                endif
              elseif @mpwr==round(@mpwr) && abs(@mpwr)<32
                zjk = (r=sqrt(r)) + flip(zj)
                if @mpwr<0.0
                  zri = r/zri
                  zjk = 1.0/zjk
                else
                  zri = zri/r
                endif
                if (abs(@mpwr)%4)>1
                  zp2 = sqr(zri)
                  z1p2 = sqr(zjk)
                endif
                if (abs(@mpwr)%8)>3
                  if (abs(@mpwr)%4)>1
                    zp4 = sqr(zp2)
                    z1p4 = sqr(z1p2)
                  else
                    zp4 = zri^4
                    z1p4 = zjk^4
                  endif
                endif
                if (abs(@mpwr)%16)>7
                  if (abs(@mpwr)%8)>3
                    zp8 = sqr(zp4)
                    z1p8 = sqr(z1p4)
                  elseif (abs(@mpwr)%4)>1
                    zp8 = zp2^4
                    z1p8 = z1p2^4
                  else
                    zp8 = sqr(zri^4)
                    z1p8 = sqr(zjk^4)
                  endif
                endif
                if abs(@mpwr)>15
                  if (abs(@mpwr)%16)>7
                    zp16 = sqr(zp8)
                    z1p16 = sqr(z1p8)
                  elseif (abs(@mpwr)%8)>3
                    zp16 = zp4^4
                    z1p16 = z1p4^4
                  elseif (abs(@mpwr)%4)>1
                    zp16 = sqr(zp2^4)
                    z1p16 = sqr(z1p2^4)
                  else
                    zp16 = (zri^4)^4
                    z1p16 = (zjk^4)^4
                  endif
                endif
                if abs(@mpwr)==31
                  zri = zp16*zp8*zp4*zp2*zri
                  zjk = z1p16*z1p8*z1p4*z1p2*zjk
                elseif abs(@mpwr)==30
                  zri = zp16*zp8*zp4*zp2
                  zjk = z1p16*z1p8*z1p4*z1p2
                elseif abs(@mpwr)==29
                  zri = zp16*zp8*zp4*zri
                  zjk = z1p16*z1p8*z1p4*zjk
                elseif abs(@mpwr)==28
                  zri = zp16*zp8*zp4
                  zjk = z1p16*z1p8*z1p4
                elseif abs(@mpwr)==27
                  zri = zp16*zp8*zp2*zri
                  zjk = z1p16*z1p8*z1p2*zjk
                elseif abs(@mpwr)==26
                  zri = zp16*zp8*zp2
                  zjk = z1p16*z1p8*z1p2
                elseif abs(@mpwr)==25
                  zri = zp16*zp8*zri
                  zjk = z1p16*z1p8*zjk
                elseif abs(@mpwr)==24
                  zri = zp16*zp8
                  zjk = z1p16*z1p8
                elseif abs(@mpwr)==23
                  zri = zp16*zp4*zp2*zri
                  zjk = z1p16*z1p4*z1p2*zjk
                elseif abs(@mpwr)==22
                  zri = zp16*zp4*zp2
                  zjk = z1p16*z1p4*z1p2
                elseif abs(@mpwr)==21
                  zri = zp16*zp4*zri
                  zjk = z1p16*z1p4*zjk
                elseif abs(@mpwr)==20
                  zri = zp16*zp4
                  zjk = z1p16*z1p4
                elseif abs(@mpwr)==19
                  zri = zp16*zp2*zri
                  zjk = z1p16*z1p2*zjk
                elseif abs(@mpwr)==18
                  zri = zp16*zp2
                  zjk = z1p16*z1p2
                elseif abs(@mpwr)==17
                  zri = zp16*zri
                  zjk = z1p16*zjk
                elseif abs(@mpwr)==16
                  zri = zp16
                  zjk = z1p16
                elseif abs(@mpwr)==15
                  zri = zp8*zp4*zp2*zri
                  zjk = z1p8*z1p4*z1p2*zjk
                elseif abs(@mpwr)==14
                  zri = zp8*zp4*zp2
                  zjk = z1p8*z1p4*z1p2
                elseif abs(@mpwr)==13
                  zri = zp8*zp4*zri
                  zjk = z1p8*z1p4*zjk
                elseif abs(@mpwr)==12
                  zri = zp8*zp4
                  zjk = z1p8*z1p4
                elseif abs(@mpwr)==11
                  zri = zp8*zp2*zri
                  zjk = z1p8*z1p2*zjk
                elseif abs(@mpwr)==10
                  zri = zp8*zp2
                  zjk = z1p8*z1p2
                elseif abs(@mpwr)==9
                  zri = zp8*zri
                  zjk = z1p8*zjk
                elseif abs(@mpwr)==8
                  zri = zp8
                  zjk = z1p8
                elseif abs(@mpwr)==7
                  zri = zp4*zp2*zri
                  zjk = z1p4*z1p2*zjk
                elseif abs(@mpwr)==6
                  zri = zp4*zp2
                  zjk = z1p4*z1p2
                elseif abs(@mpwr)==5
                  zri = zp4*zri
                  zjk = z1p4*zjk
                elseif abs(@mpwr)==4
                  zri = zp4
                  zjk = z1p4
                elseif abs(@mpwr)==3
                  zri = zp2*zri
                  zjk = z1p2*zjk
                elseif abs(@mpwr)==2
                  zri = zp2
                  zjk = z1p2
                endif
                zri = real(zjk)*zri + cri
                if @plus
                  zj = cj + imag(zjk)
                else
                  zj = cj - imag(zjk)
                endif

Well, in the next step you are going to calculate sqrt(drx*drx+dry*dry+drz*drz), so the sign of R-Ri is of no importance..

Anyway, I tried your DE in UF, and it works, unfortunately whereas the image I posted earlier took 34 seconds, the one below took more than 12 minutes, so while it is a lot more accurate ( I see the spikes now also), it is a lot less efficient. I"m trying to find a way to optimise it..

I have just finishing rendering z^2+c Mandelbulb in 8000x8000 resolution using delta DE algorithm. Step = 0.1*DE was used for very high quality. Image was downscaled to 4000x4000 for 2xAntialiasing

(http://th02.deviantart.net/fs70/300W/i/2010/010/5/f/Mandelbrot_3D___4000x4000_by_KrzysztofMarczak.jpg)
http://krzysztofmarczak.deviantart.com/art/Mandelbrot-3D-4000x4000-149897511

Wow, that's stunning Buddhi. I can even see some real fractal detail there too in certain places.

oak or maple leaf mandelbulb ;)

I have sucessfully implemented Buddhi's DE in my UltraFractal formula. As Jos observed I get the little spike things sticking out. My default DE method (the Hart, Sandin and Kauffman method) gives a very similar (but not a good, imo ;)) image. The render time was 1.58 min for a 640x640 image run on a DualCore E4500. I am still doing some optimization and testing, after which I will upload the modified formula to the UltraFractal database.