|
 |
This is my algorithmically final mandel_julia_pattern.
The posted performance graph tells it all.
It includes several innovations over the previous:
1.)Mandelbrot and Julia sets unified into a single subroutine.
2.)mandel_julia_pattern allows artistic innovation of smooth
interpolation between pure Mandelbrot and pure Julia sets
via new parameter, MJFract.
3.)My previous awful interior loops have been eliminated. Their
into static data. Their repetitive execution is now reduced to
a simple switch/case.
4.)I have added special small exponent code. This is a debatable
advantage: For instance, it now is a pretty impressive 19% faster
than mandel_pattern for Exponent=2. On the other hand, it adds an
ugly outer layer of switch/case. If you deem the latter too awful,
the default of the outer switch will handle all Exponents with no
worse than a factor of 2 slow-down. The best solution would be to
figure out how to move the outer switch inside the col loop without
harming performance too much. But a factor of 2 is a factor of 2 and
19% is 19% and these sorts of speedups are quite fragile. A better
coder than myself could probably do it.
5.)Exponents up to 255 are now treated. Of course these big exponents
take longer than smaller ones. But their presence does not handicap
the smaller ones in any way in a uniprocessor situation. If you want
to discuss problems of multi-core parallelization, I think I can add
views based on experience.
6.)I have added an exit based on the initial value of MinAbsZ2.
I do not have a scientific basis for this, but it appears to
improve logical coherence.
speed-optimized algorithm, not a code. In fact, I have edited it for clarity
so that there is some possibility of its not even compiling.
I would be glad to include this algorithm into an object oriented version
of POV.
==========================================================================
#define NEXT_SQUARE {p[2]=p[0]*p[0]-p[1]*p[1]; p[3]=2.0*p[0]*p[1]; p+=2;}
#define ACCUM_POWER(N) {DBL work=ReZToE*p[Nth1[N]]-ImZToE*p[Nth1[N]+1];
ImZToE=ReZToE*p[Nth1[N]+1]+ImZToE*p[Nth1[N]]; ReZToE=work;}
static DBL mandel_julia_pattern(double EPoint[3], int Exponent, DBL MJFract)
{
static int LastExponent=-1;
int col;
/////////////////////////////////////////////////////////////////////////////
/// The following is done only once whenever the value of Exponent changes//
/// This means that Exponent should be in a class (perhaps in TPATTERN,/////
/// perhaps in a new class that embraces both mandel_patterns and///////////
/// julia_patterns). What here are "static ints" should be members of//////
/// that class, set whenever Exponent is changed!///////////////////////////
static int NBits,NOnes; ///
static int Nth1[8]; //N.B.: Loc (in reals) within complex array! ///
if(Exponent!=LastExponent) ///
{//Exponent changed, parse out bits. ///
assert((1<Exponent)&&(Exponent<256)); ///
LastExponent = Exponent; ///
NOnes=0; ///
for(NBits=0; NBits<8; NBits++) ///
if((Exponent>>NBits)==0) ///
break; ///
else if(Exponent&(1<<NBits)) ///
Nth1[NOnes++] = 2*NBits; //2 = Stride (in DBLS) of complex ///
assert(NBits>1); ///
assert(NOnes>0); ///
} ///
/// End of hacked up version of Exponent parsing////////////////////////////
/////////////////////////////////////////////////////////////////////////////
DBL ReZ, ImZ, ReX, ImX, MinAbsZ2, ReZToE, ImZToE;
DBL ReC=0.353, ImC=0.288; //default values, un-comment 2 lines below...
ReZ = EPoint[X]; // Think of ReZ+iImZ as z
ImZ = EPoint[Y]; // Think of ReX+iImX as x, x=z in Mandelbrot
/////////////////////////////////////////////////////////////////////////////
/// Interpolate between pure Mandelbrot and pure Julia sets://///////////////
/// ///
/// MJFract = 0 Pure Mandelbrot ///
/// MJFract = 1 Pure Julia ///
/// ///
// ReC = TPat->Vals.Fractal.Coord[U]; //Uncomment this!!! ///
// ImC = TPat->Vals.Fractal.Coord[V]; //Uncomment this!!! ///
ReX = (1.0-MJFract)*ReZ + MJFract*ReC; ///
ImX = (1.0-MJFract)*ImZ + MJFract*ImC; ///
/////////////////////////////////////////////////////////////////////////////
MinAbsZ2 = ReZ*ReZ + ImZ*ImZ;
if(MinAbsZ2 > 4.0)
return(fractal_exterior_color(col, ReZ, ImZ));
DBL ZTo2ToN[16]; // ZTo2ToN[2*n]+%i*ZTo2ToN[2*n+1] will become (z^(2^n))
switch(Exponent)
{
case 2:
{
for(col = 0; col < it_max; col++)
{
double *p=ZTo2ToN;
p[0] = ReZ; // z^1
p[1] = ImZ;
NEXT_SQUARE; // z^2
ReZ = ZTo2ToN[2] + ReX;
ImZ = ZTo2ToN[3] + ImX;
DBL AbsZ2 = ReZ*ReZ + ImZ*ImZ;
if(AbsZ2 > 4.0)
return(fractal_exterior_color(col, ReZ, ImZ));
//Look, Ma, no ifs!
MinAbsZ2 = 0.5*(MinAbsZ2+AbsZ2-fabs(MinAbsZ2-AbsZ2));
}
break;
}
case 3:
{
for(col = 0; col < it_max; col++)
{
double *p=ZTo2ToN;
p[0] = ReZ; // z^1
p[1] = ImZ;
NEXT_SQUARE; // z^2
ReZ = ZTo2ToN[0]*ZTo2ToN[2] - ZTo2ToN[1]*ZTo2ToN[3] + ReX;
ImZ = ZTo2ToN[0]*ZTo2ToN[3] + ZTo2ToN[1]*ZTo2ToN[2] + ImX;
DBL AbsZ2 = ReZ*ReZ + ImZ*ImZ;
if(AbsZ2 > 4.0)
return(fractal_exterior_color(col,ReZ,ImZ));
MinAbsZ2 = 0.5*(MinAbsZ2+AbsZ2-fabs(MinAbsZ2-AbsZ2));
}
break;
}
case 4:
{
for(col = 0; col < it_max; col++)
{
double *p=ZTo2ToN;
p[0] = ReZ; // z^1
p[1] = ImZ;
NEXT_SQUARE; // z^2
NEXT_SQUARE; // z^4
ReZ = ZTo2ToN[4] + ReX;
ImZ = ZTo2ToN[5] + ImX;
DBL AbsZ2 = ReZ*ReZ + ImZ*ImZ;
if(AbsZ2 > 4.0)
return(fractal_exterior_color(col, ReZ, ImZ));
MinAbsZ2 = 0.5*(MinAbsZ2+AbsZ2-fabs(MinAbsZ2-AbsZ2));
}
break;
}
default:
for(col = 0; col < it_max; col++)
{
double *p=ZTo2ToN;
p[0] = ReZ; // z^1
p[1] = ImZ;
switch (NBits)
{//Fills up NBits complex entries of ZTo2ToN
case 8: NEXT_SQUARE; //Uneeded if Exp < 128
case 7: NEXT_SQUARE; //Uneeded if Exp < 64
case 6: NEXT_SQUARE; //Uneeded if Exp < 32
case 5: NEXT_SQUARE; //Uneeded if Exp < 16
case 4: NEXT_SQUARE; //Uneeded if Exp < 8
case 3: NEXT_SQUARE; //Uneeded if Exp < 4
default: NEXT_SQUARE; //Always needed because Exp>1
}
p=ZTo2ToN;
ReZToE = ZTo2ToN[Nth1[0] ]; //One power comes in for free!
ImZToE = ZTo2ToN[Nth1[0]+1];
switch (NOnes)
{//Fills up NBits complex entries of ZTo2ToN
case 8: ACCUM_POWER(7); //Uneeded if NOnes < 8
case 7: ACCUM_POWER(6); //Uneeded if NOnes < 7
case 6: ACCUM_POWER(5); //Uneeded if NOnes < 6
case 5: ACCUM_POWER(4); //Uneeded if NOnes < 5
case 4: ACCUM_POWER(3); //Uneeded if NOnes < 4
case 3: ACCUM_POWER(2); //Uneeded if NOnes < 3
case 2: ACCUM_POWER(1); //Uneeded if NOnes < 2
}
ReZ = ReZToE + ReX; // z <= z^Exp + x
ImZ = ImZToE + ImX;
DBL AbsZ2 = ReZ*ReZ + ImZ*ImZ;
if(AbsZ2 > 4.0)
return(fractal_exterior_color(col, ReZ, ImZ));
MinAbsZ2 = 0.5*(MinAbsZ2+AbsZ2-fabs(MinAbsZ2-AbsZ2));
}
}
return(fractal_interior_color(col, ReZ, ImZ, MinAbsZ2));
}
Post a reply to this message
|
|
