PROPOSED NEW VERSION OF mandelx_pattern

The basis of this improvement is the following algorithm.  This algorithm could
not possibly be original to me, though I did not get it from any source.

PROPOSED ALGORITHM
How to compute powers in logarithmic time.  The POV docs inaccurately imply that
computing the pth power takes p-1 multiplications.  In fact it always takes
fewer than 2*log2(p) multiplications.  Here is how that speedup works:

/////////////////////////////////////////////////////////////////
// Let z=a+ib, want zToP = z^p
//    for example, take p=37
//
// Compute z, z^2, z^4, z^8, ... by successive squaring.
// Call this sequence zTo2ToN.
//
// Write these below binary expansion of p:
//
//           1      0      0      1      0      1
//           z^32   z^16   z^8    z^4    z^2    z^1
//    zToP = z^32                *z^4          *z^1  = z^37
//
// Notice that this takes:
//    5 complex multiplies to get the sequence zTo2ToN
//    2 complex multiplies accumulating -> zNew
//  ----------------------
//    7 complex multiplies total
//
// i,e., WAY less than the 36 required if computed via
//
//           z^1, z^2, z^3,..., z^34, z^35
//
// In general, to compute a pth power, it takes
//
//    CountZeros(p)+2*CountOnes(p)-2 complex multiplies
//
// For certain values of p (like 15 or 30), there are methods
// that are somewhat faster, but the extra code burden would
// probably not be worth it.
/////////////////////////////////////////////////////////////////



================================================================================
This obviates a LOT of code.....

I propose that the code below replaces:
    1.) mandelx_pattern
    2.) mandel_pattern
    3.) mandel3_pattern
    4.) mandel4_pattern
    5.) MANDEL_PATTERN
    6.) MANDEL3_PATTERN
    7.) MANDEL4_PATTERN

and all that machinery.

The same trick applies to the Julia Set as well, NEW_juliax_pattern eliminating:
    8.) juliax_pattern
    9.) julia_pattern
    10.) julia3_pattern
    11.) julia4_pattern
    12.) JULIA_PATTERN
    13.) JULIA3_PATTERN
    14.) JULIA4_PATTERN

and all that machinery.

When these are both done, there is no further need for:
    15.) binomial_coeff
    16.) InitializeBinomialCoefficients
    17.) BinomialCoefficients
    18.) BinomialCoefficientsInited

and all that machinery.

================================================================================
And it goes a LOT faster.  The enclosed figure, ALGO_TIMES.JPG, shows the
speedups.
Since I cannot post it to this forum, I will place it in some more permissive
forum on the POVRAY site.

================================================================================
Finally, it is a LOT more accurate.  This is mostly a matter of peace-of-mind
for coders.  I doubt very much that users will notice any artistic improvement
from this, but, for exponent=30, there actually were two pixels different even
on my 401x401 hacked up grid (but no smaller exponent differed at all).


================================================================================
================================================================================
Here is a hacked-up version of the code.  I have never operated it in situ (for
some reason [missing DLL?] I am unable to execute my POV compilations).  But I
am pretty sure that it works.

/********NEW_mandelx_pattern begins********************************************
#define CMPLX_SQ(A,B)      {DBL w = A*A - B*B; B = 2.0*A*B;   A = w;}
#define CMPLX_MUL(A,B,C,D) {DBL w = A*C - B*D; B = A*D + B*C; A = w;}

static DBL NEW_mandelx_pattern(double EPoint[3])
{//  assert(exponent>0);
  int col,ExpShifted;
  DBL ReZ, ImZ, ReX, ImX, MinAbsZ2;

  ReZ = ReX = EPoint[X];  // Think of ReZ+iImZ as z
  ImZ = ImX = EPoint[Y];  // Think of ReX+iImX as x, start with z=x
  MinAbsZ2 = ReZ*ReZ + ImZ*ImZ;

  for(col = 0; col < it_max; col++)
  {
      DBL ReZTo2ToN = ReZ; // Will become z^(2^bit_pos)
      DBL ImZTo2ToN = ImZ;

      for(ExpShifted=exponent; (ExpShifted&1)==0; ExpShifted >>= 1)
          CMPLX_SQ(ReZTo2ToN,ImZTo2ToN); // Next trailing zero in exponent

      DBL ReZToE = ReZTo2ToN; // Least significant 1 in exponent,
      DBL ImZToE = ImZTo2ToN; // This will become z^exponent

      for(ExpShifted >>= 1   ; ExpShifted!=0    ; ExpShifted >>= 1)
      {// Next bit-position in exponent
          CMPLX_SQ(ReZTo2ToN,ImZTo2ToN); // Square previous z^(2^bit_pos)
          if((ExpShifted&1)!=0)          // This bit-position = 1?
              CMPLX_MUL(ReZToE,ImZToE,ReZTo2ToN,ImZTo2ToN); // -> z^exponent
      }

      ReZ = ReZToE + ReX; // z <= z^exponent + x
      ImZ = ImZToE + ImX;

      DBL AbsZ2 = ReZ*ReZ + ImZ*ImZ;

      if(AbsZ2 < MinAbsZ2) MinAbsZ2 = AbsZ2;
      if(AbsZ2 > 4.0)
          return(fractal_exterior_color(col, ReZ, ImZ));
  }

  return(fractal_interior_color(col, ReZ, ImZ, MinAbsZ2));
}
********NEW_mandelx_pattern ends**********************************************/

Post a reply to this message