Welcome to Fractal Forums

z[0, _, _] := 0
z[n_, c_, -1] := {1/2 (1 - Sqrt[1 - 4 z[n - 1, c, -1]]) + c, 
                  1/2 (1 + Sqrt[1 - 4 z[n - 1, c, -1]]) + c}
z[n_, c_,  1] := {1/2 (1 - Sqrt[1 - 4 z[n - 1, c,  1]]) + c, 
                  1/2 (1 + Sqrt[1 - 4 z[n - 1, c,  1]]) + c}
z[n_, c_] := {z[n, c, -1], z[n, c, 1]}

import collections
from cmath import sqrt


def flatten(l):
    for el in l:
        if isinstance(el, collections.Iterable) and not isinstance(el, (str, bytes)):
            yield from flatten(el)
        else:
            yield el


def z(n, c, b=None):
    if n == 0:
        yield 0
    elif b == -1:
        for _z in flatten(z(n-1, c, -1)):
            yield .5 * (1 - sqrt(1 - 4 * _z)) + c
            yield .5 * (1 + sqrt(1 - 4 * _z)) + c
    elif b == 1:
        for _z in flatten(z(n - 1, c, 1)):
            yield .5 * (1 - sqrt(1 - 4 * _z)) + c
            yield .5 * (1 + sqrt(1 - 4 * _z)) + c
    elif b is None:
        yield flatten(z(n, c, -1))
        yield flatten(z(n, c, 1))

Ian1-M {
;
; Ian McDonald custom formula from the forum (Mandelbrot mode)
;
init:
  z = (0,0)
  c = #pixel
loop:
  z = ((((z*z+c)*(z*z+c)+(c*c+z))*((z*z+c)*(z*z+c)+(c*c+z))+((c*c+z)*(c*c+z)+(z*z+c)))*(((z*z+c)*(z*z+c)+(c*c+z))*((z*z+c)*(z*z+c)+(c*c+z))+((c*c+z)*(c*c+z)+(z*z+c)))+(((c*c+z)*(c*c+z)+(z*z+c))*((c*c+z)*(c*c+z)+(z*z+c))+((z*z+c)*(z*z+c)+(c*c+z))))
bailout:
  |z| <= @bailout
default:
  title = "Ian1 Mandelbrot"
  float param bailout
    caption = "Bailout value"
    default = 4.0
  endparam
switch:
  type = "Ian1-J"
  seed = #pixel
  bailout = bailout
}

Ian1-J {
;
; Ian McDonald custom formula from the forum (Julia mode)
;
init:
  z = #pixel
  c = @seed
loop:
  z = ((((z*z+c)*(z*z+c)+(c*c+z))*((z*z+c)*(z*z+c)+(c*c+z))+((c*c+z)*(c*c+z)+(z*z+c)))*(((z*z+c)*(z*z+c)+(c*c+z))*((z*z+c)*(z*z+c)+(c*c+z))+((c*c+z)*(c*c+z)+(z*z+c)))+(((c*c+z)*(c*c+z)+(z*z+c))*((c*c+z)*(c*c+z)+(z*z+c))+((z*z+c)*(z*z+c)+(c*c+z))))
bailout:
  |z| <= @bailout
default:
  title = "Ian1 Julia"
  param seed
    caption = "Seed"
    default = (-0.1, 0)
  endparam
  float param bailout
    caption = "Bailout value"
    default = 4.0
  endparam
switch:
  type = "Ian1-M"
  bailout = bailout
}

for (i = 0; i < Iterations; i++) {
 z = (
   cMul(
	(
	    cMul((cMul((cMul(z,z)+c),(cMul(z,z)+c))+(cMul(c,c)+z)),
		 (cMul((cMul(z,z)+c),(cMul(z,z)+c))+(cMul(c,c)+z)))+
			(cMul((cMul(c,c)+z),(cMul(c,c)+z))+(cMul(z,z)+c))
	),
	(
	    cMul((cMul((cMul(z,z)+c),(cMul(z,z)+c))+(cMul(c,c)+z)),
		 (cMul((cMul(z,z)+c),(cMul(z,z)+c))+(cMul(c,c)+z)))+
			(cMul((cMul(c,c)+z),(cMul(c,c)+z))+(cMul(z,z)+c))
	)
	)+
	(
	    cMul((cMul((cMul(c,c)+z),(cMul(c,c)+z))+(cMul(z,z)+c)),
		 (cMul((cMul(c,c)+z),(cMul(c,c)+z))+(cMul(z,z)+c)))+
			(cMul((cMul(z,z)+c),(cMul(z,z)+c))+(cMul(c,c)+z))
	)
    );
		if (dot(z,z)> 100.0) break;
}

Power = 4;

vec2 Plot(vec2 a) {
  vec2 b;
  b.x = sqrt (a.x * a.x + a.y * a.y);
  b.y = atan (a.y / a.x);
  if ( a.x < 0.) b.y += 3.141592;
  b.x = pow(b.x,float(Power));
  b.y *= float(Power);
  a.x = cos(b.y) * b.x;
  a.y = sin(b.y) * b.x;
  return a;
}

...
      for (i = 0; i < Iterations; i++) {
        z += Plot(z) + c;
        if (dot(z,z)> 100.0) break;
      }
...

import random as rnd
import sys as sys
import threading as thrd
from cmath import sqrt
from math import floor
from PIL import Image
import numpy as np


class Tree:
    def __init__(self, v0, fl, fr):
        self.v = v0
        self.fl = fl
        self.fr = fr

    def left(self):
        return Tree(self.fl(self.v), self.fl, self.fr)

    def right(self):
        return Tree(self.fr(self.v), self.fl, self.fr)


def empty():  # generator without elements.
    return    # looks like a dirty hack.
    yield     # as far as I know the only way.


def chain(*its):
    for it in its:
        for i_ in it:
            yield i_


def bf(t):  # traverses infinite tree breadth first such that it gets to all the nodes and gives all the values eventually. For n iterations, need 2**(n+1)-1 points.
    ts = (t for t in [t])  # generator with one element
    while ts:
        nts = empty()
        for t in ts:
            if t is not None:
                yield t.v
                nts = chain(nts, [t.left()], [t.right()])
        ts = nts


  # I decided to sample from a circular region. Could go with rectangular region or with something more carefully chosen to get more detail in the busy central area

def rangle():  # generate random numbers evenly distributed on a circle in the complex plane
    while True:
        z = complex(rnd.uniform(-1, 1), rnd.uniform(-1, 1))
        yield z*z/(z*z.conjugate())


def rcircle(c0, r0): # generate random numbers evenly distributed in a circle in the complex plane
    a = rangle()
    while True:
        r = sqrt(rnd.uniform(0, r0*r0))
        yield r * next(a) - c0


def mtree(c_, z0_):  # corecursively generates all the points to be plotted
    return Tree(z0_, lambda z: .5*(1-sqrt(1-4*z))+c_, lambda z: .5*(1+sqrt(1-4*z))+c_)


def mbrot(iterations, resolution, filename='MBrot'):

    filename += '.png'

    # triple code for the three colors

    histr = np.array([[0. for _ in range(resolution)] for _ in range(resolution)], dtype=np.float32)
    histg = np.array([[0. for _ in range(resolution)] for _ in range(resolution)], dtype=np.float32)
    histb = np.array([[0. for _ in range(resolution)] for _ in range(resolution)], dtype=np.float32)
    img = Image.fromarray(np.uint8(histr)*255)
    img.save(filename)

    cs = rcircle(.0+.0j, 2.)

    itr = 2**(iterations + 1) - 1  # perhaps this shouldn't be hard-coded. As is, green will have half the number of points as red and blue will have half that
    itg = 2**iterations - 1
    itb = 2**(iterations - 1) - 1

    try:
        while True:
            tree = mtree(next(cs), next(cs))  # I decided to go Buddhagram-style, also giving a random initial z0. Could fix. Perhaps at the critical point z = .25+0i.
            points = bf(tree)

            for z, it_ in zip(points, range(itr)):
                r = (z.real + 3.)/6.  # the +3/6 determines what point is in the middle of the image (as is the 0,0 point) and how large the overall area is. Perhaps rewrite in terms of center point and zoom level and don't hard-code.
                x = floor((1. - 1.e-10) * resolution * r)  # the  1.e-10 is to avoid being right on the edge of a pixel
                if 0 <= x < resolution:
                    im = (z.imag + 3.)/6.
                    y = floor((1. - 1.e-10) * resolution * im)
                    if 0 <= y < resolution:
                        histr[x, y] += 1
                        if it_ <= itg:
                            histg[x, y] += 1
                            if it_ <= itb:
                                histb[x, y] += 1

                        im = (3. - z.imag)/6
                    y = floor((1.-1.e-10) * resolution * im)  # double code for enforced mirror-symmetry
                    if 0 <= y < resolution:
                        histr[x, y] += 1
                        if it_ <= itg:
                            histg[x, y] += 1
                            if it_ <= itb:
                                histb[x, y] += 1
            imhistr, binsr = np.histogram(histr.flatten(), 2 ** 16, normed=True)  # I'm not even sure if the 2 ** 16 even makes a difference
            imhistg, binsg = np.histogram(histr.flatten(), 2 ** 16, normed=True) # perhaps 2 ** 8 = 256 suffices.
            imhistb, binsb = np.histogram(histr.flatten(), 2 ** 16, normed=True) # It's how finely the histogram should be divided for equalization.
            cdfr = imhistr.cumsum()
            cdfr = cdfr / cdfr[-1]
            imeqr = np.interp(histr.flatten(), binsr[:-1], cdfr)
            mineqr = np.min(imeqr)
            maxeqr = np.max(imeqr)
            if mineqr != maxeqr:
                imeqr -= mineqr
                imeqr *= 255/(maxeqr-mineqr)
            cdfg = imhistg.cumsum()
            cdfg = cdfg / cdfg[-1]
            imeqg = np.interp(histg.flatten(), binsg[:-1], cdfg)
            mineqg = np.min(imeqg)
            maxeqg = np.max(imeqg)
            if mineqg != maxeqg:
                imeqg -= mineqg
                imeqg *= 255/(maxeqg-mineqg)
            cdfb = imhistb.cumsum()
            cdfb = cdfb / cdfb[-1]
            imeqb = np.interp(histb.flatten(), binsb[:-1], cdfb)
            mineqb = np.min(imeqb)
            maxeqb = np.max(imeqb)
            if mineqb != maxeqb:
                imeqb -= mineqb
                imeqb *= 255/(maxeqb-mineqb)
            imgr = Image.fromarray(imeqr.reshape(histr.shape)).convert('L')  # perhaps these four lines could be made into a single step 
            imgg = Image.fromarray(imeqg.reshape(histg.shape)).convert('L')
            imgb = Image.fromarray(imeqb.reshape(histb.shape)).convert('L')
            img = Image.merge('RGB', [imgr, imgg, imgb])
            try:  # if I don't do this, it crashes when something else accesses the file.
                img.save(filename)
            except IOError:
                pass
    except KeyboardInterrupt:  # break loop when pressing Ctrl+C
        img.save(filename)  # save a final time

def runbrot():  # need a parameter-less wrapper function to allow for deeper nesting -> higher iterations
    it = int(input("Please enter an iteration depth.
"
                   "(Note: If this is too large, it will crash. Try something in the very low double-digits.
"
                   "It will also be exponentially slower."))
    res = int(input("Please Enter a resolution."))
    fn = input("How would you like to call your file? (Without file ending)")
    print("To end, press Ctrl+C")
    mbrot(it, res, fn)

sys.setrecursionlimit(2 ** 28)  # lift Python-interpreter-imposed nesting limit
thrd.stack_size(2**28-1)  # lift OS-imposed nesting limit - will crash when out of memory!
thread = thrd.Thread(target=runbrot)
thread.start()

#version 400

uniform vec2 resolution;
uniform vec3 palette[256];
uniform double xmin;
uniform double xmax;
uniform double ymin;
uniform double ymax;
uniform double bailout;
uniform int maxiters;
uniform int samplepixels;

double sqrsamplepixels=double(samplepixels*samplepixels);
double bailout_squared=double(bailout*bailout);
double magnitude,r1,r2,g1,g2,b1,b2,tweenval,x,y;
double realiters;
vec4 finalcol,col;
int superx,supery;
double stepx=(xmax-xmin)/resolution.x/double(samplepixels);
double stepy=(ymax-ymin)/resolution.y/double(samplepixels);
int colval,colval2;
dvec2 z,c,lastz,lastz2,s1,s2,n,d,a,d2,temp,temp2,temp3;

dvec2 cmul( dvec2 a, dvec2 b )  { return dvec2( a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x ); }
dvec2 csqr( dvec2 a ) { return dvec2(a.x*a.x-a.y*a.y, 2.0*a.x*a.y ); }

void main(void)
{
	finalcol=vec4(0,0,0,0);

	//C is constant
	c.x = -0.1;
	c.y = 0.0;

	for (supery=0;supery<samplepixels;supery++)
	{
		for (superx=0;superx<samplepixels;superx++)
		{
			//z from current pixel
			z.x = xmin+gl_FragCoord.x/resolution.x*(xmax-xmin)+(stepx*double(superx));
			z.y = ymin+gl_FragCoord.y/resolution.y*(ymax-ymin)+(stepy*double(supery));
			int i;

			for(i=0; i<maxiters; i++) 
			{
				// http://www.fractalforums.com/images-showcase-(rate-my-fractal)/meta-mandelbrot-(mandeljulia)/?PHPSESSID=b5fd43c3acfd4c3472a92524915f0100
				dvec2 ztmp=cmul((cmul((csqr((csqr(z)+c))+(csqr(c)+z)),(csqr((csqr(z)+c))+(csqr(c)+z)))+(csqr((csqr(c)+z))+(csqr(z)+c))),(cmul((csqr((csqr(z)+c))+(csqr(c)+z)),(csqr((csqr(z)+c))+(csqr(c)+z)))+(csqr((csqr(c)+z))+(csqr(z)+c))))+(cmul((csqr((csqr(c)+z))+(csqr(z)+c)),(csqr((csqr(c)+z))+(csqr(z)+c)))+(csqr((csqr(z)+c))+(csqr(c)+z)));

				magnitude=(ztmp.x*ztmp.x+ztmp.y*ztmp.y);
				if (magnitude>bailout_squared) break;
				z=ztmp;
			}

			if (i==maxiters) {
				col=vec4(0.0,0.0,0.0,1.0);
			} else {
				//CPM smooth colors
				//note that double precision does not support log so it needs to be cast as float
				realiters=float(i+1-((log(log(sqrt(float(magnitude))))/log(16)))); //normally log is the power factor, needs to be raised here to smooth colors
				colval=int(mod(realiters,255));
				colval2=int(mod(colval+1,255));
				tweenval=realiters-int(realiters);
				r1=palette[colval].r;
				g1=palette[colval].g;
				b1=palette[colval].b;
				r2=palette[colval2].r;
				g2=palette[colval2].g;
				b2=palette[colval2].b;
				col=vec4(r1+((r2-r1)*tweenval),g1+((g2-g1)*tweenval),b1+((b2-b1)*tweenval),1.0);
			}

			finalcol+=col;
		}
	}
	gl_FragColor = vec4(finalcol/double(sqrsamplepixels));
}
	

#define LEVEL0(a,b) (a)
#define LEVEL1(a,b) csqr(LEVEL0(a,b))+LEVEL0(b,a)
#define LEVEL2(a,b) csqr(LEVEL1(a,b))+LEVEL1(b,a)
#define LEVEL3(a,b) csqr(LEVEL2(a,b))+LEVEL2(b,a)
#define LEVEL4(a,b) csqr(LEVEL3(a,b))+LEVEL3(b,a)
dvec2 ztmp = LEVEL4(z,c);

(((((z^P1)+(z^P3+c))*((z^P1)+(z^P3+c))+((c^P2)+(c^P4+z)))*(((z^P1)+(z^P3+c))*((z^P1)+(z^P3+c))+((c^P2)+(c^P4+z)))+(((c^P2)+(c^P4+z))*((c^P2)+(c^P4+z))+((z^P1)+(z^P3+c))))*((((z^P1)+(z^P3+c))*((z^P1)+(z^P3+c))+((c^P2)+(c^P4+z)))*(((z^P1)+(z^P3+c))*((z^P1)+(z^P3+c))+((c^P2)+(c^P4+z)))+(((c^P2)+(c^P4+z))*((c^P2)+(c^P4+z))+((z^P1)+(z^P3+c))))+((((c^P2)+(c^P4+z))*((c^P2)+(c^P4+z))+((z^P1)+(z^P3+c)))*(((c^P2)+(c^P4+z))*((c^P2)+(c^P4+z))+((z^P1)+(z^P3+c)))+(((z^P1)+(z^P3+c))*((z^P1)+(z^P3+c))+((c^P2)+(c^P4+z)))))