Mandelbrotai {
;
; Mandelbrot that uses a modified (branch-free) log: re(logz)+i*abs(im(logz))
; The downside: lots of streching, even on integer powers.
init:
  $define DEBUG
  z = @start
  complex rep = real(@power);
  complex iimp = 1i*imag(@power);
loop:
  ;z = z^@power + #pixel
  if z==0
       z=#pixel;
  else
      complex logz = log(z);
      complex logaz = real(logz)+1i*abs(imag(logz));
      ;complex logbz = log(real(z))+1i*abs(imag(z)); even wierder.
      z = exp(rep*logaz+iimp*logaz)+#pixel;
  endif
bailout:
  |z| <= @bailout
default:
  title = "Special log mandelbrot"
  center = (-0.5, 0)
$IFDEF VER50
  rating = recommended
$ENDIF
  param start
    caption = "Starting point"
    default = (0,0)
    hint = "The starting point parameter can be used to distort the Mandelbrot \
            set. Use (0, 0) for the standard MandelbrotAI set."
  endparam
  param power
    caption = "Power"
    default = (2,0)
    hint = "This parameter sets the exponent for the Mandelbrot formula. \
            Increasing the real part to 3, 4, and so on, will add stalks to \
            the Mandelbrot figure. You are always stuck with distortion. Use (2, 0) \
            for the standard MandelbrotAI set."
  endparam
  float param bailout
    caption = "Bailout value"
    default = 1e8
    min = 1.0
$IFDEF VER40
    exponential = true
$ENDIF
    hint = "This parameter defines how soon an orbit bails out while \
            iterating. Larger values give smoother outlines; values around 4 \
            give more interesting shapes around the set. Values much less than 4 \
            will distort the fractal."
  endparam
switch:
  type = "MandelbrotaiJulia"
  seed = #pixel
  power = power
  bailout = bailout
}

MandelbrotaiJulia {
;
; Julia version of Special log mandelbrot.
;
init:
  $define DEBUG
  z = #pixel
  complex rep = real(@power);
  complex iimp = 1i*imag(@power);
loop:
  ;z = z^@power + #pixel
  if z==0
       z=@seed;
  else
      complex logz = log(z);
      complex logaz = real(logz)+1i*abs(imag(logz));
      ;complex logbz = log(real(z))+1i*abs(imag(z)); even wierder.
      z = exp(rep*logaz+iimp*logaz)+@seed;
  endif
bailout:
  |z| <= @bailout
default:
  title = "Special log Julia"
  center = (-0.5, 0)
$IFDEF VER50
  rating = recommended
$ENDIF
  param power
    caption = "Power"
    default = (2,0)
    hint = "This parameter sets the exponent for the Julia formula. \
            Increasing the real part to 3, 4, and so on, will add discs to \
            the Julia figure. No distortion on the Julia. Use (2, 0) \
            for the standard JuliaAI set."
  endparam
  param seed
    caption = "Julia seed"
    default = (-1.25, 0)
    hint = "Use this parameter to create many different Julia sets. A good \
            way to set this parameter is with the Switch, Eyedropper, or \
            Explore features."
  endparam
  float param bailout
    caption = "Bailout value"
    default = 1e8
    min = 1.0
$IFDEF VER40
    exponential = true
$ENDIF
    hint = "This parameter defines how soon an orbit bails out while \
            iterating. Larger values give smoother outlines; values around 4 \
            give more interesting shapes around the set. Values much less than 4 \
            will distort the fractal."
  endparam
}

rootChoose {
;Idea independantly thought of as: http://www.righto.com/2011/12/new-multi-branch-algorithm-to-render.html
;Square roots are a problem: there are two choices. We choose both!
;generalized mandelbrot set, (z=z^(A+B/N)+c, A,B,N are integers) with this rule:
;at EVERY iteration, we branch out to all possible choices of root (N branches)
;The combinatorial explosion isn't that bad near the perephery because
;if you don't pick the exact right combination you will go flying off in most
;places and the tree is pruned quickly.
;What about complex exponents? They have infinite #'s of branches, and some cause escape to be impossible.
      ;fortunatly, we can invert the imaginary component and get two that glue up.
      ;this causes distortion, so we lack a better solution.
;However, in the denser regions it is a problem, increase "max branch" if you
;see broken lines.
init:
     ;$define DEBUG

     ;TODO: history. complex history[#maxiter+1]; best history.
     float leaveFraction = 0; fraction staying in the branch.

     int nBranchLeft = @maxBranchedIter; set a limit on how many we can do.

     complex zBranch[#maxIter+1];
     ;setLength(zBranch,#maxIter);
     zBranch[0] = @start; #pixel;

     int choiceBranch[#maxIter];which root we choose at [ind] to determine zBranch[ind+1].
     ;setLength(zBranch,#maxIter);
     int ind=0;
     while ind<#maxIter-1; fill with zeros (this IS nessessary?)
           choiceBranch[ind] = 0;
           ind = ind+1;
     endwhile



     int denomP = @denomPower;
     int numerP = @numPower;
     int intP = @intPower;
     float denomF = @denomPower;  the float version (avoids int-int errors).
     float bailout2 = @bailout;
     bool calcPorp = @calcPorp;
     int placeValue=0; which z-value we are on.
     float porpInside=0;how much we discovered is inside.
     float porpOutside=0;how much we discovered is outside.
     int maxIter=0;the maximum iterations we do.
     complex maxIterZ=0;
     float imagP = @imagPower;
     int IBM=1; stands for imaginary branch multiplier.
     if imagP !=0,
        IBM = IBM*2;
     endif
     while nBranchLeft>0; DONT USE A RECURSIVE FUNCTION! RECURSION = MEMORY CRASH!
           bool backUp=false; wheter we diverged and thus backed out of the tree.
           complex z = zBranch[placeValue];
           float weight = denomF^(-placeValue);
           if placeValue==#maxIter, ;no further to go!
              backUp=true;
              porpInside = porpInside+weight;
              if calcPorp==false, ;dont bother continuing.
                 nBranchLeft = -1;
              endif
           elseif |z|>bailout2; flung away.
              backUp=true;
              porpOutside = porpOutside+weight;
           endif
           if backUp==true,;BACK UP: take the next branch, if any. Otherwise set nBranchLeft <0 to quit.
              int leaf = placeValue-1;
              if leaf==-1, ;nessessary (only?) for points out of bail to begin with.
                  nBranchLeft = -10;BREAK out of the loop (and why no break statement in UF5?).
              else; pick the next branch.
                  choiceBranch[leaf] = choiceBranch[leaf]+1;change the branch that got us here.
                  while choiceBranch[leaf]==denomP*IBM && nBranchLeft>0;carry the numbers to the larger value.
                        choiceBranch[leaf]=0; roll over to 0.
                        choiceBranch[leaf-1] = choiceBranch[leaf-1]+1; carry.
                        leaf=leaf-1;
                      if choiceBranch[0]==denomP*IBM, ;we ran out of branches.
                          nBranchLeft= -1;
                      endif
                  endwhile
              endif
              placeValue = leaf+1;
           elseif nBranchLeft>0; don't back up: continue on our current path.
              placeValue = placeValue+1;
              if placeValue>maxIter,;record set. placeValue is bounded to maxIter.
                   maxIter = placeValue;
                   maxIterZ = z;
              endif
           endif

           ;we UPDATED place value, now z doesn't make sense.
           ;instead we use the parent from one place above.
           if nBranchLeft>0 && placeValue<#maxIter, ;if nothing in the backUp loop caused us to exit.
              ;recalculate z at placeValue based on z before placeValue and the formula from
              ;http://www.mathamazement.com/Lessons/Pre-Calculus/06_Additional-Topics-in-Trigonometry/de-moivres-theorem.html.
                  ;Lemma 6.6.4
              complex zParent = zBranch[placeValue-1];
              int option = choiceBranch[placeValue-1]; which one we choose, from 0 to IBM*denomP-1.
              int roption=option;
              if option>=denomP,;the real value of the option (ignore the second branch on the logrythim).
                   roption=option-denomP;
              endif
              int ioption=0;imaginary component.
              if option>=denomP,
                 ioption=1;
              endif
              complex newZ=0;
              if denomP==2 && (numerP+2*intP==4);optimization for 2+0/2
                 newZ = zParent*zParent*(1-2*roption)+#pixel;
              elseif denomP==2 && (numerP+2*intP==5);optimization for 2+1/2
                 newZ = zParent*zParent*sqrt(zParent)*(1-2*roption);
              else
                   float cR = cabs(zParent);
                   float cTheta = 0; convert into polar.
                   if cR>0, ;if r is zero, theta doesn't matter, but avoid the singularity.
                      cTheta = atan2(zParent);
                   endif
                   float newR = cR^(intP+numerP/denomP);
                   float newTheta = cTheta*(intP+numerP/denomP)+(2*#pi*option)/denomP;
                   newZ=newR*cos(newTheta)+newR*1i*sin(newTheta);
              endif

              ;add in the imaginary exponent:
              if imagP !=0 && |newZ|>0,
                  complex logParent = log(zParent);
                  if ioption==1,;branch by multiplying the imaginary by -1.
                       logParent = real(logParent)-1i*imag(logParent);
                  endif
                  complex fac = exp(1i*imagP*logParent);a^b = exp(b*log(a))
                  newZ = newZ*fac;
                  ;complex logaz = real(log(newZ))+1i*abs(imag(log(newZ)));
                  ;newZ = newZ*exp(logaz*imagP*1i);special absolute log value.
              endif
              zBranch[placeValue] = newZ+#pixel;
           endif
           nBranchLeft = nBranchLeft-1; limit total sub-branches.
     endwhile

     if porpInside==0 && porpOutside==0,;avoid indeterminate singularity.
        porpInside = 1;
        porpOutside = 1;
     endif
    float stayFraction = porpInside/(porpInside+porpOutside);
loop: ;just an empty loop to the same number of steps as above.
    maxIter=maxIter-1;
    if stayFraction>0, ;interior coloring.
       z = stayFraction;
    else ;exterior coloring.
       z=maxIterZ;
    endif
bailout:
    maxIter>=0
default:
  title = "Which square root?"
  center = (0, 0)
  int param intPower
    caption = "Integer part of exponent";
    hint = "this is for convience, since the numerator can also be changed, i.e. 1+1/2 = 0+3/2"
    default = 2
  endparam
  int param numPower
    caption = "Numerator exponant"
    hint = "Try zero!"
    default = 1
  endparam
  int param denomPower
    caption = "Denominator exponant"
    hint = "This is the number of branch points, no matter what the other real values are."
    default = 2
  endparam
  float param imagPower
    caption = "Imaginary exponant"
    hint = "Uses the absolute value of the log's imag to avoid branches, but introduces distortion"
    default = 0
  endparam
  bool param calcPorp
    caption = "calculate fraction inside"
    hint = "Turn on if you want to use the cooresponding inside coloring algorythim (but it is slower)."
    default = false
  endparam
  float param bailout
    caption = "Bailout"
    hint = "The standard bailout, <3.5 or so = cut off fractal.";
    default = 16
    $IFDEF VER40; what does this do?
           exponential = true
    $ENDIF
  endparam
  int param maxBranchedIter
     caption = "max branch";
     default=20000;
  endparam
  maxiter = 200;
  ;doesn't seem to work: if #calculationPurpose==1
  ;     maxiter=12;
  ;endif
  periodicity = 0; causes problems!
  param start
    caption = "Starting point"
    default = (0,0)
    hint = "The starting point parameter can be used to distort the Mandelbrot \
            set. Use (0, 0) for the standard MandelbrotAI set."
  endparam
switch:
  type = "rootChooseJulia"
  seed = #pixel
  intPower = intPower
  numPower = numPower
  denomPower = denomPower
  imagPower = imagPower
  maxBranchedIter = maxBranchedIter
  bailout = bailout
}

rootChooseJulia{
;Julia version of "Which square root?"
init:
     ;$define DEBUG

     z = #pixel;

     ;TODO: history. complex history[#maxiter+1]; best history.
     float leaveFraction = 0; fraction staying in the branch.

     int nBranchLeft = @maxBranchedIter; set a limit on how many we can do.

     complex zBranch[#maxIter+1];
     ;setLength(zBranch,#maxIter);
     zBranch[0] = #pixel;

     int choiceBranch[#maxIter];which root we choose at [ind] to determine zBranch[ind+1].
     ;setLength(zBranch,#maxIter);
     int ind=0;
     while ind<#maxIter-1; fill with zeros (this IS nessessary?)
           choiceBranch[ind] = 0;
           ind = ind+1;
     endwhile

     int denomP = @denomPower;
     int numerP = @numPower;
     int intP = @intPower;
     float denomF = @denomPower;  the float version (avoids int-int errors).
     float bailout2 = @bailout;
     bool calcPorp = @calcPorp;
     int placeValue=0; which z-value we are on.
     float porpInside=0;how much we discovered is inside.
     float porpOutside=0;how much we discovered is outside.
     int maxIter=0;the maximum iterations we do.
     complex maxIterZ=0;
     float imagP = @imagPower;
     int IBM=1; stands for imaginary branch multiplier.
     if imagP !=0,
        IBM = IBM*2;
     endif
     while nBranchLeft>0; DONT USE A RECURSIVE FUNCTION! RECURSION = MEMORY CRASH!
           bool backUp=false; wheter we diverged and thus backed out of the tree.
           complex z = zBranch[placeValue];
           float weight = denomF^(-placeValue);
           if placeValue==#maxIter, ;no further to go!
              backUp=true;
              porpInside = porpInside+weight;
              if calcPorp==false, ;dont bother continuing.
                 nBranchLeft = -1;
              endif
           elseif |z|>bailout2; flung away.
              backUp=true;
              porpOutside = porpOutside+weight;
           endif
           if backUp==true,;BACK UP: take the next branch, if any. Otherwise set nBranchLeft <0 to quit.
              int leaf = placeValue-1;
              if leaf==-1, ;nessessary (only?) for points out of bail to begin with.
                  nBranchLeft = -10;BREAK out of the loop (and why no break statement in UF5?).
              else; pick the next branch.
                  choiceBranch[leaf] = choiceBranch[leaf]+1;change the branch that got us here.
                  while choiceBranch[leaf]==denomP*IBM && nBranchLeft>0;carry the numbers to the larger value.
                        choiceBranch[leaf]=0; roll over to 0.
                        choiceBranch[leaf-1] = choiceBranch[leaf-1]+1; carry.
                        leaf=leaf-1;
                      if choiceBranch[0]==denomP*IBM, ;we ran out of branches.
                          nBranchLeft= -1;
                      endif
                  endwhile
              endif
              placeValue = leaf+1;
           elseif nBranchLeft>0; don't back up: continue on our current path.
              placeValue = placeValue+1;
              if placeValue>maxIter,;record set. placeValue is bounded to maxIter.
                   maxIter = placeValue;
                   maxIterZ = z;
              endif
           endif

           ;we UPDATED place value, now z doesn't make sense.
           ;instead we use the parent from one place above.
           if nBranchLeft>0 && placeValue<#maxIter, ;if nothing in the backUp loop caused us to exit.
              ;recalculate z at placeValue based on z before placeValue and the formula from
              ;http://www.mathamazement.com/Lessons/Pre-Calculus/06_Additional-Topics-in-Trigonometry/de-moivres-theorem.html.
                  ;Lemma 6.6.4
              complex zParent = zBranch[placeValue-1];
              int option = choiceBranch[placeValue-1]; which one we choose, from 0 to IBM*denomP-1.
              int roption=option;
              if option>=denomP,;the real value of the option (ignore the second branch on the logrythim).
                   roption=option-denomP;
              endif
              int ioption=0;imaginary component.
              if option>=denomP,
                 ioption=1;
              endif
              complex newZ=0;
              if denomP==2 && (numerP+2*intP==4);optimization for 2+0/2
                 newZ = zParent*zParent*(1-2*roption)+@seed;
              elseif denomP==2 && (numerP+2*intP==5);optimization for 2+1/2
                 newZ = zParent*zParent*sqrt(zParent)*(1-2*roption);
              else
                   float cR = cabs(zParent);
                   float cTheta = 0; convert into polar.
                   if cR>0, ;if r is zero, theta doesn't matter, but avoid the singularity.
                      cTheta = atan2(zParent);
                   endif
                   float newR = cR^(intP+numerP/denomP);
                   float newTheta = cTheta*(intP+numerP/denomP)+(2*#pi*option)/denomP;
                   newZ=newR*cos(newTheta)+newR*1i*sin(newTheta);
              endif

              ;add in the imaginary exponent:
              if imagP !=0 && |newZ|>0,
                  complex logParent = log(zParent);
                  if ioption==1,;branch by multiplying the imaginary by -1.
                       logParent = real(logParent)-1i*imag(logParent);
                  endif
                  complex fac = exp(1i*imagP*logParent);a^b = exp(b*log(a))
                  newZ = newZ*fac;
                  ;complex logaz = real(log(newZ))+1i*abs(imag(log(newZ)));
                  ;newZ = newZ*exp(logaz*imagP*1i);special absolute log value.
              endif
              zBranch[placeValue] = newZ+@seed;
           endif
           nBranchLeft = nBranchLeft-1; limit total sub-branches.
     endwhile

     if porpInside==0 && porpOutside==0,;avoid indeterminate singularity.
        porpInside = 1;
        porpOutside = 1;
     endif
    float stayFraction = porpInside/(porpInside+porpOutside);
loop: ;just an empty loop to the same number of steps as above.
    maxIter=maxIter-1;
    if stayFraction>0, ;interior coloring.
       z = stayFraction;
    else ;exterior coloring.
       z=maxIterZ;
    endif
bailout:
    maxIter>=0
default:
  title = "Which square root (Julia)?"
  center = (0, 0)
  int param intPower
    caption = "Integer part of exponent";
    hint = "this is for convience, since the numerator can also be changed, i.e. 1+1/2 = 0+3/2"
    default = 2
  endparam
  int param numPower
    caption = "Numerator exponant"
    hint = "Try zero!"
    default = 1
  endparam
  int param denomPower
    caption = "Denominator exponant"
    hint = "This is the number of branch points, no matter what the other real values are."
    default = 2
  endparam
  float param imagPower
    caption = "Imaginary exponant"
    hint = "Uses the absolute value of the log's imag to avoid branches, but introduces distortion"
    default = 0
  endparam
  bool param calcPorp
    caption = "calculate fraction inside"
    hint = "Turn on if you want to use the cooresponding inside coloring algorythim (but it is slower)."
    default = false
  endparam
  param seed
    caption = "Julia seed"
    default = (-1.25, 0)
    hint = "Use this parameter to create many different Julia sets. A good \
            way to set this parameter is with the Switch, Eyedropper, or \
            Explore features."
  endparam
  float param bailout
    caption = "Bailout"
    hint = "The standard bailout";
    default = 16
    $IFDEF VER40; what does this do?
           exponential = true
    $ENDIF
  endparam
  int param maxBranchedIter
     caption = "max branch";
     default=20000;
  endparam
  maxiter = 200;
  ;doesn't seem to work: if #calculationPurpose==1
  ;     maxiter=12;
  ;endif
  periodicity = 0; causes problems!
}