Fractal Software > Programming

Complex functions implementation by 2 reals

(1/1)

Alef:
This could be usefull for writing formulas. I found this usefull when I needed complex power 3 by 2 real variables. Wikipedia have some formulas but in more powers you 'll need to do algebra by yourself and that is quite cumbersome.

--- Quote ---Formula Compiler tutorial
Introduction to the Fractal Explorer Formula Compiler
version 2.03.02
September 2004
Written by Sirotinsky Arthur, Olga Fedorenko
and Denis McCauley

--- End quote ---

4.7 The complex functions

You can use any function of the complex variables in your formula, but you must describe them, like this (remember, that Z = X + jY and C = CR + jCI):

z*z:
X_new := (X-Y)*(X+Y);
Y_new := 2 * X*Y;

z*c:
X_new := X*CR - Y*CI;
Y_new := X*CI + Y*CR;

z*z*z:
X_new := X*(X*X-3*Y*Y);                       // -this is optimized formula
Y_new := Y*(3*X*X-Y*Y);

z*z*z*z:
tmp_r := (X-Y)*(X+Y);
tmp_i := 2 * X*Y;
X_new := (tmp_r - tmp_i) * (tmp_r + tmp_i);
Y_new := 2 * tmp_r * tmp_i;

1/z:
tmp   := X*X + Y*Y + 1E-25;                   // -without 1E-25 may by
X_new := X / tmp;                             //  divizion by zero
Y_new :=-Y / tmp;

1/(z*z):
tmp   := X*X + Y*Y;
tmp   := tmp * tmp + 1E-25;                   // -see above
X_new := (X*X - Y*Y) / tmp;
Y_new := (2*X * Y  ) / tmp;

Sqrt(z):
tmp   := Sqrt(X*X + Y*Y);
X_new := Sqrt(Abs((X + tmp)/2));
Y_new := Sqrt(Abs((tmp - X)/2));

Exp(z):
tmp   := Exp(X);                              // e^z
X_new := tmp * Cos(Y);
Y_new := tmp * Sin(Y);

Exp(1/z):
tmp   := X*X + Y*Y + 1E-25;                   // -see above
s1    := X/tmp;
s2    :=-Y/tmp;
tmp   := Exp(s1);
X_new := tmp * Cos(s2);
Y_new := tmp * Sin(s2);

Ln(z):
X_new :=Log2(X*X+Y*Y)/2.7182818285;
Y_new :=ArcTan2(Y,X);

z^c:
h1x   :=Log2(X*X+Y*Y)/2.7182818285;           // -Ln(z)
h1y   :=ArcTan2(Y,X);
h2x   :=h1x*CR - h1y*CI;                      // -Ln(z)*c
h2y   :=h1y*CR + h1x*CI;
f     :=Exp(h2x);                             // -Exp(Ln(z)*c)
X_new :=f*Cos(h2y);
Y_new :=f*Sin(h2y);

Sin(z):
tmp   := Exp(Y)/2;                            // -optimized formula (!)
tmp2  := 0.25/tmp;
X_new := Sin(X) * (tmp+tmp2);
Y_new := Cos(X) * (tmp-tmp2);

Cos(z):
X_new := Cos(X)*Cosh(Y);                      // -not optimized formula
Y_new :=-Sin(X)*Sinh(Y);

Tan(z):
X_new := Sin(2*X)/(Cos(2*X)+Cosh(2*Y));
Y_new := Sin(2*Y)/(Cos(2*X)+Cosh(2*Y));

Sinh(z):                      // -hiperbolic sinus
X_new := Sinh(X)*Cos(Y);                      // -it is not optimized formula
Y_new := Cosh(X)*Sin(Y);

Cosh(z):                      // -hiperbolic cosinus
X_new := Cosh(X)*Cos(Y);                      // -it is not optimized formula
Y_new := Sinh(X)*Sin(Y);

claude:
for correct results, strongly suggest NOT adding a constant to avoid division by 0.  no constant will be small enough for deep zooms. instead check for 0 explicitly and/or normalize NaNs/infinities afterwards (e.g. with C99 cproj() function).  note: correct behaviour will probably will be slower

hobold:
I am wondering in which cases (of 1/z) the addition of a tiny epsilon makes an actual difference. I mean when is the tiny constant not effectively rounded to zero, because the left value (sum of squares) is large enough that the right addend is underflowing during alignment?

Hmm, then again, the sum of squares itself could be tiny, and then the supposedly tiny constant may be larger ... destroying any information that was in the input numbers.

In other words, one would need to either add a much smaller constant, which in turn could still cause the division to overflow, or explicitly check the magnitude before deciding on a suitable value for a constant. In the latter case, one can check for zero in the first place and handle the problematic case separately.

Hmm yeah ... there seems to be no simple solution that is behaved better over a larger range of values.

Alef:
I think it have to do with compiler and alsou initial values. I guess floating variable never can realy get value zero by calculation. But I don't have theoretical knowledge about this;) Alsou I had not done deep zooms, maybe with petrurbation it eventualy could introduce distortion. But it could help to owercome compilers error messages;)

In Chaos Pro programm shows messages like "Fractal (Z*Z+C) may be the result of arithmetic exceptions. The calculation function produced invalid results for at least 7 pixels." Most probably it generated infinities.

xenodreambuie:
I think the question of how to handle divide by zero depends on several factors, such as how deep you want to go and by what methods, what precision you are using, what is supported/convenient in your programming environment, and also what the users expect. Most users don't care how many pixels may have been technically NaN somewhere, they just want a nice picture.

In any case, adding a constant is not a good idea. There are other problems with those simple methods of invert and divide: intermediate precision loss due to overflow and underflow. With better methods you need to test anyway, so it costs little more to explicitly test for 0. Here is a method for invert to minimize precision loss, in Pascal.

--- Code: ---function CInvert(const c1:TComplex):TComplex;
var R,Denom: Extended;
begin
if Abs(c1.x)>=Abs(c1.y) then begin
if c1.x<>0 then begin
R:= c1.y/c1.x;
Denom:= 1/(c1.x+R*c1.y);
Result.x:= Denom;
Result.y:= - R*Result.x;
end else begin
result.x:= 1e100; // or handle it in some other way
Result.y:= -1e100;
end;
end else begin
R:= c1.x/c1.y;
Denom:= -1/(c1.y+R*c1.x);
Result.y:= Denom;
Result.x:= -R*result.y;
end;
end;

--- End code ---
If you're just using real x and y instead of a complex type, it's easily enough changed.