import enum
import time
import numba as nb
import numpy as np
import matplotlib
from PIL import Image, PyAccess
# Amount of times to print the total progress
PROGRESS_STEPS: int = 20
# Number of pixels (width*height) and aspect ratio (width/height)
RESOLUTION: int = 1920*1080
ASPECT_RATIO: float = 1920/1080
# Value of the center pixel
CENTER: complex = 0 + 0j # For testing: -4.6875 + 2.63671875j
# Value range of the real part (width of the horizontal axis)
RE_RANGE: float = 10 # For testing: 10/16
# Show grid lines for integer real and imaginary parts
SHOW_GRID: bool = False
GRID_COLOR: tuple[int, int, int] = (125, 125, 125)
# Matplotlib named colormap
COLORMAP_NAME: str = 'inferno'
# Color of the interior of the fractal (convergent points)
INTERIOR_COLOR: tuple[int, int, int] = (0, 0, 60)
# Color for large divergence counts
ERROR_COLOR: tuple[int, int, int] = (0, 0, 60)
# Plot range of the axes
X_MIN = CENTER.real - RE_RANGE/2 # min Re(z)
X_MAX = CENTER.real + RE_RANGE/2 # max Re(z)
Y_MIN = CENTER.imag - RE_RANGE/(2*ASPECT_RATIO) # min Im(z)
Y_MAX = CENTER.imag + RE_RANGE/(2*ASPECT_RATIO) # max Im(z)
x_range = X_MAX - X_MIN
y_range = Y_MAX - Y_MIN
pixels_per_unit = np.sqrt(RESOLUTION/(x_range*y_range))
# Width and height of the image in pixels
WIDTH = round(pixels_per_unit*x_range)
HEIGHT = round(pixels_per_unit*y_range)
# Maximum iterations for the divergence test
MAX_ITER: int = 10**4 # recommended >= 10**3
# Minimum consecutive abs(r) decreases to declare linear divergence
MIN_R_DROPS: int = 4 # recommended >= 2
# Minimum iterations to start checking for slow drift (unknown divergence)
MIN_ITER_SLOW: int = 200 # recommended >= 100
# Max value of Re(z) and Im(z) such that the recursion doesn't overflow
CUTOFF_RE = 7.564545572282618e+153
CUTOFF_IM = 112.10398935569289
# Precompute the colormap
CMAP_LEN: int = 2000
cmap_mpl = matplotlib.colormaps[COLORMAP_NAME]
# Start away from 0 (discard black values for the 'inferno' colormap)
# Matplotlib's colormaps have 256 discrete color points
n_cmap = round(256*0.85)
CMAP = [cmap_mpl(k/256) for k in range(256 - n_cmap, 256)]
# Interpolate
x = np.linspace(0, 1, num=CMAP_LEN)
xp = np.linspace(0, 1, num=n_cmap)
c0, c1, c2 = tuple(np.interp(x, xp, [c[k] for c in CMAP]) for k in range(3))
CMAP = []
for x0, x1, x2 in zip(c0, c1, c2):
CMAP.append(tuple(round(255*x) for x in (x0, x1, x2)))
class DivType(enum.Enum):
"""Divergence type."""
MAX_ITER = 0 # Maximum iterations reached
SLOW = 1 # Detected slow growth (maximum iterations will be reached)
CYCLE = 2 # Cycled back to the same value after 8 iterations
LINEAR = 3 # Detected linear divergence
CUTOFF_RE = 4 # Diverged by exceeding the real part cutoff
CUTOFF_IM = 5 # Diverged by exceeding the imaginary part cutoff
@nb.jit(nb.float64(nb.float64, nb.int64), nopython=True)
def smooth(x, k=1):
"""Recursive exponential smoothing function."""
y = np.expm1(np.pi*x)/np.expm1(np.pi)
if k <= 1:
return y
return smooth(y, np.fmin(6, k - 1))
@nb.jit(nb.float64(nb.float64), nopython=True)
def get_delta_im(x):
"""Get the fractional part of the smoothed divergence count for
imaginary part blow-up."""
nu = np.log(np.abs(x)/CUTOFF_IM)/(np.pi*CUTOFF_IM - np.log(CUTOFF_IM))
nu = np.fmax(0, np.fmin(nu, 1))
return smooth(1 - nu, 2)
@nb.jit(nb.float64(nb.float64, nb.float64), nopython=True)
def get_delta_re(x, e):
"""Get the fractional part of the smoothed divergence count for
real part blow-up."""
nu = np.log(np.abs(x)/CUTOFF_RE)/np.log1p(e)
nu = np.fmax(0, np.fmin(nu, 1))
return 1 - nu
@nb.jit(
nb.types.containers.Tuple((
nb.float64,
nb.types.EnumMember(DivType, nb.int64)
))(nb.complex128),
nopython=True
)
def divergence_count(z):
"""Return a smoothed divergence count and the type of divergence."""
delta_im = -1
delta_re = -1
cycle = 0
r0 = -1
r_drops = 0 # Counter for number of consecutive times abs(r) decreases
a, b = z.real, z.imag
a_cycle, b_cycle = a, b
cutoff_re_squared = CUTOFF_RE*CUTOFF_RE
for k in range(MAX_ITER):
e = 0.5*np.exp(-np.pi*b)
cycle += 1
if cycle == 8:
cycle = 0
r = e*np.hypot(0.5 + a, b)/(1e-6 + np.abs(b))
if r < r0 < 0.5:
r_drops += 1
else:
r_drops = 0
# Stop early due to likely slow linear divergence
if r_drops >= MIN_R_DROPS:
delta = 0.25*(CUTOFF_RE - a)
return k + delta, DivType.LINEAR
# Detected slow growth (maximum iterations will be reached)
if ((k >= MIN_ITER_SLOW) and (r0 <= r)
and (r + (r - r0)*(MAX_ITER - k) < 8*0.05)):
delta = 0.25*(CUTOFF_RE - a)
return k + delta, DivType.SLOW
r0 = r
# Cycled back to the same value after 8 iterations
if (a - a_cycle)**2 + (b - b_cycle)**2 < 1e-16:
return k, DivType.CYCLE
a_cycle = a
b_cycle = b
a0 = a
# b0 = b
s = np.sin(np.pi*a)
c = np.cos(np.pi*a)
# Equivalent to:
# z = 0.25 + z - (0.25 + 0.5*z)*np.exp(np.pi*z*1j)
# where z = a + b*1j
a += e*(b*s - (0.5 + a)*c) + 0.25
b -= e*(b*c + (0.5 + a0)*s)
if b < -CUTOFF_IM:
delta_im = get_delta_im(-b)
if a*a + b*b > cutoff_re_squared:
delta_re = get_delta_re(np.hypot(a, b), e)
# Diverged by exceeding a cutoff
if delta_im >= 0 or delta_re >= 0:
if delta_re < 0 or delta_im <= delta_re:
return k + delta_im, DivType.CUTOFF_IM
else:
return k + delta_re, DivType.CUTOFF_RE
# Maximum iterations reached
return -1, DivType.MAX_ITER
@nb.jit(nb.complex128(nb.float64, nb.float64), nopython=True)
def pixel_to_z(a, b):
re = X_MIN + (X_MAX - X_MIN)*a/WIDTH
im = Y_MAX - (Y_MAX - Y_MIN)*b/HEIGHT
return re + 1j*im
@nb.jit(nb.float64(nb.float64), nopython=True)
def cyclic_map(g):
"""A continuous function that cycles back and forth between 0 and 1."""
# This can be any continuous function.
# Log scale removes high-frequency color cycles.
# freq_div = 1
# g = np.log1p(np.fmax(0, g/freq_div)) - np.log1p(1/freq_div)
# Beyond this value for float64, decimals are truncated
if g >= 2**51:
return -1
# Normalize and cycle
# g += 0.5 # phase from 0 to 1
return 1 - np.abs(2*(g - np.floor(g)) - 1)
def get_pixel(px, py):
z = pixel_to_z(px, py)
dc, div_type = divergence_count(z)
match div_type:
case DivType.MAX_ITER
