Data update
This commit is contained in:
Ingy döt Net
parent
4d5544505c
commit
4924dd0264
3073 changed files with 55820 additions and 4408 deletions
|
|
@ -1,5 +1,123 @@
|
|||
from functools import reduce
|
||||
import numba
|
||||
numba.config.CUDA_ENABLE_PYNVJITLINK = True # prevent cuda ptx version errors
|
||||
|
||||
def mandelbrot(x, y, c): return ' *'[abs(reduce(lambda z, _: z*z + c, range(99), 0)) < 2]
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
print('\n'.join(''.join(mandelbrot(x, y, x/50 + y/50j) for x in range(-100, 25)) for y in range(-50, 50)))
|
||||
import cupy as cp
|
||||
import numba.cuda as cuda
|
||||
|
||||
import decimal as dc # decimal floating point arithmetic with arbitrary precision
|
||||
dc.getcontext().prec = 80 # set precision to 80 digits (about 256 bits)
|
||||
|
||||
d, h = 1600, 1000 # pixel density (= image width) and image height
|
||||
n, r = 100000, 100000.0 # number of iterations and escape radius (r > 2)
|
||||
|
||||
a = dc.Decimal("-1.256827152259138864846434197797294538253477389787308085590211144291")
|
||||
b = dc.Decimal(".37933802890364143684096784819544060002129071484943239316486643285025")
|
||||
|
||||
S = np.zeros(n+1, dtype=np.complex128)
|
||||
u, v = dc.Decimal(0), dc.Decimal(0)
|
||||
|
||||
for i in range(n+1):
|
||||
S[i] = float(u) + float(v) * 1j
|
||||
if u * u + v * v < r * r:
|
||||
u, v = u * u - v * v + a, 2 * u * v + b
|
||||
else:
|
||||
print("The reference sequence diverges within %s iterations." % i)
|
||||
break
|
||||
|
||||
x = np.linspace(0, 2, num=d+1, dtype=np.float64)
|
||||
y = np.linspace(0, 2 * h / d, num=h+1, dtype=np.float64)
|
||||
|
||||
A, B = np.meshgrid(x - 1, y - h / d)
|
||||
C = 5.0e-35 * (A + B * 1j)
|
||||
|
||||
def iteration_cupy_cuda(S, C):
|
||||
I = cp.zeros(C.shape, dtype=np.int32)
|
||||
E, Z, dZ = cp.zeros_like(C), cp.zeros_like(C), cp.zeros_like(C)
|
||||
|
||||
iteration = cp.RawKernel("""
|
||||
#include <cupy/complex.cuh>
|
||||
|
||||
extern "C" __global__
|
||||
void iterate(int dim_x, int dim_y, int n, double r,
|
||||
complex<double> *S, complex<double> *C,
|
||||
int *I, complex<double> *E, complex<double> *Z, complex<double> *dZ) {
|
||||
|
||||
int x = blockIdx.x * blockDim.x + threadIdx.x;
|
||||
int y = blockIdx.y * blockDim.y + threadIdx.y;
|
||||
|
||||
if (x < dim_x and y < dim_y) { // prevent memory access errors
|
||||
int x_y = x * dim_y + y; // cupy arrays are in row-major order
|
||||
complex<double> delta = C[x_y];
|
||||
|
||||
int index = I[x_y];
|
||||
complex<double> epsilon = E[x_y];
|
||||
complex<double> z = Z[x_y];
|
||||
complex<double> dz = dZ[x_y];
|
||||
|
||||
double abs2_r = r * r;
|
||||
double abs2_z, abs2_e;
|
||||
|
||||
for (int i = 0; i < n; i++) {
|
||||
abs2_z = z.real() * z.real() + z.imag() * z.imag();
|
||||
abs2_e = epsilon.real() * epsilon.real() + epsilon.imag() * epsilon.imag();
|
||||
|
||||
if (abs2_z < abs2_e) { // rebase when z is closer to zero
|
||||
epsilon = z; index = 0; // reset reference orbit
|
||||
}
|
||||
if (abs2_z < abs2_r) {
|
||||
epsilon = (2. * S[index] + epsilon) * epsilon + delta; index = index + 1;
|
||||
dz = 2. * z * dz + 1.; z = S[index] + epsilon;
|
||||
}
|
||||
else {
|
||||
break;
|
||||
}
|
||||
}
|
||||
I[x_y] = index; E[x_y] = epsilon; Z[x_y] = z; dZ[x_y] = dz;
|
||||
}
|
||||
}
|
||||
""", "iterate")
|
||||
|
||||
griddim, blockdim = (h // 32 + 1, d // 32 + 1), (32, 32)
|
||||
iteration(griddim, blockdim, (h+1, d+1, n, r, cp.asarray(S), cp.asarray(C), I, E, Z, dZ))
|
||||
return I.get(), E.get(), Z.get(), dZ.get()
|
||||
|
||||
def iteration_numba_cuda(S, C):
|
||||
I = cp.zeros(C.shape, dtype=np.int32)
|
||||
E, Z, dZ = cp.zeros_like(C), cp.zeros_like(C), cp.zeros_like(C)
|
||||
|
||||
@cuda.jit()
|
||||
def iteration(S, C, I, E, Z, dZ):
|
||||
x, y = cuda.grid(2)
|
||||
|
||||
if x < h+1 and y < d+1: # prevent memory access errors
|
||||
delta, index, epsilon, z, dz = C[x, y], I[x, y], E[x, y], Z[x, y], dZ[x, y]
|
||||
|
||||
def abs2(z):
|
||||
return z.real * z.real + z.imag * z.imag
|
||||
|
||||
for i in range(n):
|
||||
if abs2(z) < abs2(epsilon): # rebase when z is closer to zero
|
||||
index, epsilon = 0, z # reset reference orbit
|
||||
if abs2(z) < abs2(r):
|
||||
index, epsilon = index + 1, (2 * S[index] + epsilon) * epsilon + delta
|
||||
z, dz = S[index] + epsilon, 2 * z * dz + 1
|
||||
else:
|
||||
break
|
||||
|
||||
I[x, y], E[x, y], Z[x, y], dZ[x, y] = index, epsilon, z, dz
|
||||
|
||||
griddim, blockdim = (h // 32 + 1, d // 32 + 1), (32, 32)
|
||||
iteration[griddim, blockdim](cp.asarray(S), cp.asarray(C), I, E, Z, dZ)
|
||||
return I.get(), E.get(), Z.get(), dZ.get()
|
||||
|
||||
I, E, Z, dZ = iteration_numba_cuda(S, C) # use iteration_cupy_cuda or iteration_numba_cuda
|
||||
D = np.zeros(C.shape, dtype=np.float64)
|
||||
|
||||
N = abs(Z) > 2 # exterior distance estimation
|
||||
D[N] = np.log(abs(Z[N])) * abs(Z[N]) / abs(dZ[N])
|
||||
|
||||
plt.imshow(D ** 0.15, cmap=plt.cm.jet, origin="lower")
|
||||
plt.savefig("Mandelbrot_deep_zoom.png", dpi=300)
|
||||
|
|
|
|||
8
Task/Mandelbrot-set/Python/mandelbrot-set-11.py
Normal file
8
Task/Mandelbrot-set/Python/mandelbrot-set-11.py
Normal file
|
|
@ -0,0 +1,8 @@
|
|||
print(
|
||||
'\n'.join(
|
||||
''.join(
|
||||
' *'[(z:=0, c:=x/50+y/50j, [z:=z*z+c for _ in range(99)], abs(z))[-1]<2]
|
||||
for x in range(-100,25)
|
||||
)
|
||||
for y in range(-50,50)
|
||||
))
|
||||
5
Task/Mandelbrot-set/Python/mandelbrot-set-12.py
Normal file
5
Task/Mandelbrot-set/Python/mandelbrot-set-12.py
Normal file
|
|
@ -0,0 +1,5 @@
|
|||
from functools import reduce
|
||||
|
||||
def mandelbrot(x, y, c): return ' *'[abs(reduce(lambda z, _: z*z + c, range(99), 0)) < 2]
|
||||
|
||||
print('\n'.join(''.join(mandelbrot(x, y, x/50 + y/50j) for x in range(-100, 25)) for y in range(-50, 50)))
|
||||
|
|
@ -10,13 +10,23 @@ y = np.linspace(0, 2 * h / d, num=h+1)
|
|||
A, B = np.meshgrid(x - 1, y - h / d)
|
||||
C = 2.0 * (A + B * 1j) - 0.5
|
||||
|
||||
Z, dZ = np.zeros_like(C), np.zeros_like(C)
|
||||
D, S, T = np.zeros(C.shape), np.zeros(C.shape), np.zeros(C.shape)
|
||||
def iteration(C):
|
||||
S, T = np.zeros(C.shape), np.zeros(C.shape)
|
||||
Z, dZ = np.zeros_like(C), np.zeros_like(C)
|
||||
|
||||
for k in range(n):
|
||||
M = abs(Z) < r
|
||||
S[M], T[M] = S[M] + np.exp(- abs(Z[M])), T[M] + 1
|
||||
Z[M], dZ[M] = Z[M] ** 2 + C[M], 2 * Z[M] * dZ[M] + 1
|
||||
def iterate(C, S, T, Z, dZ):
|
||||
S, T = S + np.exp(- abs(Z)), T + 1
|
||||
Z, dZ = Z * Z + C, 2 * Z * dZ + 1
|
||||
return S, T, Z, dZ
|
||||
|
||||
for i in range(n):
|
||||
M = abs(Z) < r
|
||||
S[M], T[M], Z[M], dZ[M] = iterate(C[M], S[M], T[M], Z[M], dZ[M])
|
||||
|
||||
return S, T, Z, dZ
|
||||
|
||||
S, T, Z, dZ = iteration(C)
|
||||
D = np.zeros(C.shape)
|
||||
|
||||
plt.imshow(S ** 0.1, cmap=plt.cm.twilight_shifted, origin="lower")
|
||||
plt.savefig("Mandelbrot_set_1.png", dpi=200)
|
||||
|
|
|
|||
|
|
@ -13,13 +13,23 @@ y = np.linspace(0, 2 * h / d, num=h+1)
|
|||
A, B = np.meshgrid(x - 1, y - h / d)
|
||||
C = (2.0 + 1.0j) * (A + B * 1j) - 0.5
|
||||
|
||||
Z, dZ, ddZ = np.zeros_like(C), np.zeros_like(C), np.zeros_like(C)
|
||||
D, S, T = np.zeros(C.shape), np.zeros(C.shape), np.zeros(C.shape)
|
||||
def iteration(C):
|
||||
S, T = np.zeros(C.shape), np.zeros(C.shape)
|
||||
Z, dZ, ddZ = np.zeros_like(C), np.zeros_like(C), np.zeros_like(C)
|
||||
|
||||
for k in range(n):
|
||||
M = abs(Z) < r
|
||||
S[M], T[M] = S[M] + np.sin(density * np.angle(Z[M])), T[M] + 1
|
||||
Z[M], dZ[M], ddZ[M] = Z[M] ** 2 + C[M], 2 * Z[M] * dZ[M] + 1, 2 * (dZ[M] ** 2 + Z[M] * ddZ[M])
|
||||
def iterate(C, S, T, Z, dZ, ddZ):
|
||||
S, T = S + np.sin(density * np.angle(Z)), T + 1
|
||||
Z, dZ, ddZ = Z * Z + C, 2 * Z * dZ + 1, 2 * (dZ * dZ + Z * ddZ)
|
||||
return S, T, Z, dZ, ddZ
|
||||
|
||||
for i in range(n):
|
||||
M = abs(Z) < r
|
||||
S[M], T[M], Z[M], dZ[M], ddZ[M] = iterate(C[M], S[M], T[M], Z[M], dZ[M], ddZ[M])
|
||||
|
||||
return S, T, Z, dZ, ddZ
|
||||
|
||||
S, T, Z, dZ, ddZ = iteration(C)
|
||||
D = np.zeros(C.shape)
|
||||
|
||||
N = abs(Z) >= r # basic normal map effect and stripe average coloring (potential function)
|
||||
P, Q = S[N] / T[N], (S[N] + np.sin(density * np.angle(Z[N]))) / (T[N] + 1)
|
||||
|
|
|
|||
|
|
@ -4,8 +4,8 @@ import matplotlib.pyplot as plt
|
|||
d, h = 200, 1200 # pixel density (= image width) and image height
|
||||
n, r = 8000, 10000 # number of iterations and escape radius (r > 2)
|
||||
|
||||
a = -.743643887037158704752191506114774 # https://mathr.co.uk/web/m-location-analysis.html
|
||||
b = 0.131825904205311970493132056385139 # try: a, b, n = -1.748764520194788535, 3e-13, 800
|
||||
a = -.743643887037158704752191506114774 # coordinates by github.com/josch
|
||||
b = 0.131825904205311970493132056385139 # https://github.com/josch/mandelbrot
|
||||
|
||||
x = np.linspace(0, 2, num=d+1)
|
||||
y = np.linspace(0, 2 * h / d, num=h+1)
|
||||
|
|
@ -13,12 +13,21 @@ y = np.linspace(0, 2 * h / d, num=h+1)
|
|||
A, B = np.meshgrid(x * np.pi, y * np.pi)
|
||||
C = 8.0 * np.exp((A + B * 1j) * 1j) + (a + b * 1j)
|
||||
|
||||
Z, dZ = np.zeros_like(C), np.zeros_like(C)
|
||||
D = np.zeros(C.shape)
|
||||
def iteration(C):
|
||||
Z, dZ = np.zeros_like(C), np.zeros_like(C)
|
||||
|
||||
for k in range(n):
|
||||
M = Z.real ** 2 + Z.imag ** 2 < r ** 2
|
||||
Z[M], dZ[M] = Z[M] ** 2 + C[M], 2 * Z[M] * dZ[M] + 1
|
||||
def iterate(C, Z, dZ):
|
||||
Z, dZ = Z * Z + C, 2 * Z * dZ + 1
|
||||
return Z, dZ
|
||||
|
||||
for i in range(n):
|
||||
M = abs(Z) < r
|
||||
Z[M], dZ[M] = iterate(C[M], Z[M], dZ[M])
|
||||
|
||||
return Z, dZ
|
||||
|
||||
Z, dZ = iteration(C)
|
||||
D = np.zeros(C.shape)
|
||||
|
||||
N = abs(Z) > 2 # exterior distance estimation
|
||||
D[N] = np.log(abs(Z[N])) * abs(Z[N]) / abs(dZ[N])
|
||||
|
|
@ -26,12 +35,14 @@ D[N] = np.log(abs(Z[N])) * abs(Z[N]) / abs(dZ[N])
|
|||
plt.imshow(D.T ** 0.05, cmap=plt.cm.nipy_spectral, origin="lower")
|
||||
plt.savefig("Mercator_Mandelbrot_map.png", dpi=200)
|
||||
|
||||
X, Y = C.real, C.imag # zoom images (adjust circle size 100 and zoom level 20 as needed)
|
||||
R, c, z = 100 * (2 / d) * np.pi * np.exp(- B), min(d, h) + 1, max(0, h - d) // 20
|
||||
M = 50 * (2 / d) * np.pi * np.exp(- B) # adjust marker size 50 as needed
|
||||
k, l = min(d, h) + 1, max(0, h - d) // 20 # adjust zoom level 20 as needed
|
||||
|
||||
fig, axs = plt.subplots(2, 3, figsize=(12, 8))
|
||||
for i, ax in enumerate(axs.flat):
|
||||
X, Y = C[i*l:i*l+k, 0:d].real, C[i*l:i*l+k, 0:d].imag
|
||||
S, T = M[0:k, 0:d] ** 2, D[i*l:i*l+k, 0:d] ** 0.5
|
||||
ax.scatter(X, Y, s=S, c=T, cmap=plt.cm.nipy_spectral)
|
||||
ax.axis('equal')
|
||||
|
||||
fig, ax = plt.subplots(2, 2, figsize=(12, 12))
|
||||
ax[0,0].scatter(X[1*z:1*z+c,0:d], Y[1*z:1*z+c,0:d], s=R[0:c,0:d]**2, c=D[1*z:1*z+c,0:d]**.5, cmap=plt.cm.nipy_spectral)
|
||||
ax[0,1].scatter(X[2*z:2*z+c,0:d], Y[2*z:2*z+c,0:d], s=R[0:c,0:d]**2, c=D[2*z:2*z+c,0:d]**.4, cmap=plt.cm.nipy_spectral)
|
||||
ax[1,0].scatter(X[3*z:3*z+c,0:d], Y[3*z:3*z+c,0:d], s=R[0:c,0:d]**2, c=D[3*z:3*z+c,0:d]**.3, cmap=plt.cm.nipy_spectral)
|
||||
ax[1,1].scatter(X[4*z:4*z+c,0:d], Y[4*z:4*z+c,0:d], s=R[0:c,0:d]**2, c=D[4*z:4*z+c,0:d]**.2, cmap=plt.cm.nipy_spectral)
|
||||
plt.savefig("Mercator_Mandelbrot_zoom.png", dpi=100)
|
||||
|
|
|
|||
|
|
@ -1,12 +1,13 @@
|
|||
import numba
|
||||
|
||||
import numpy as np
|
||||
import numba as nb
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
import decimal as dc # decimal floating point arithmetic with arbitrary precision
|
||||
dc.getcontext().prec = 80 # set precision to 80 digits (about 256 bits)
|
||||
|
||||
d, h = 100, 2000 # pixel density (= image width) and image height
|
||||
n, r = 80000, 100000 # number of iterations and escape radius (r > 2)
|
||||
n, r = 80000, 100000.0 # number of iterations and escape radius (r > 2)
|
||||
|
||||
a = dc.Decimal("-1.256827152259138864846434197797294538253477389787308085590211144291")
|
||||
b = dc.Decimal(".37933802890364143684096784819544060002129071484943239316486643285025")
|
||||
|
|
@ -14,12 +15,12 @@ b = dc.Decimal(".379338028903641436840967848195440600021290714849432393164866432
|
|||
S = np.zeros(n+1, dtype=np.complex128)
|
||||
u, v = dc.Decimal(0), dc.Decimal(0)
|
||||
|
||||
for k in range(n+1):
|
||||
S[k] = float(u) + float(v) * 1j
|
||||
for i in range(n+1):
|
||||
S[i] = float(u) + float(v) * 1j
|
||||
if u * u + v * v < r * r:
|
||||
u, v = u * u - v * v + a, 2 * u * v + b
|
||||
else:
|
||||
print("The reference sequence diverges within %s iterations." % k)
|
||||
print("The reference sequence diverges within %s iterations." % i)
|
||||
break
|
||||
|
||||
x = np.linspace(0, 2, num=d+1, dtype=np.float64)
|
||||
|
|
@ -28,37 +29,37 @@ y = np.linspace(0, 2 * h / d, num=h+1, dtype=np.float64)
|
|||
A, B = np.meshgrid(x * np.pi, y * np.pi)
|
||||
C = (- 8.0) * np.exp((A + B * 1j) * 1j)
|
||||
|
||||
@nb.njit(parallel=True)
|
||||
def calculation(C):
|
||||
E, I = np.zeros_like(C), np.zeros(C.shape, dtype=np.int64)
|
||||
Z, dZ = np.zeros_like(C), np.zeros_like(C)
|
||||
@numba.njit(parallel=True)
|
||||
def iteration_numba(S, C):
|
||||
I = np.zeros(C.shape, dtype=np.intp)
|
||||
E, Z, dZ = np.zeros_like(C), np.zeros_like(C), np.zeros_like(C)
|
||||
|
||||
def iteration(C):
|
||||
E, I = np.zeros_like(C), np.zeros(C.shape, dtype=np.int64)
|
||||
Z, dZ = np.zeros_like(C), np.zeros_like(C)
|
||||
def iteration(S, C):
|
||||
I = np.zeros(C.shape, dtype=np.intp)
|
||||
E, Z, dZ = np.zeros_like(C), np.zeros_like(C), np.zeros_like(C)
|
||||
|
||||
def abs2(z):
|
||||
return z.real * z.real + z.imag * z.imag
|
||||
|
||||
def iterate(E, I, Z, dZ, C):
|
||||
E, I = (2 * S[I] + E) * E + C, I + 1
|
||||
def iterate(C, I, E, Z, dZ):
|
||||
I, E = I + 1, (2 * S[I] + E) * E + C
|
||||
Z, dZ = S[I] + E, 2 * Z * dZ + 1
|
||||
return E, I, Z, dZ
|
||||
return I, E, Z, dZ
|
||||
|
||||
for k in range(n):
|
||||
M = abs2(Z) < abs2(E)
|
||||
E[M], I[M] = Z[M], 0 # rebase when z is closer to zero
|
||||
for i in range(n):
|
||||
M = abs2(Z) < abs2(E) # rebase when z is closer to zero
|
||||
I[M], E[M] = 0, Z[M] # reset the reference orbit
|
||||
M = abs2(Z) < abs2(r)
|
||||
E[M], I[M], Z[M], dZ[M] = iterate(E[M], I[M], Z[M], dZ[M], C[M])
|
||||
I[M], E[M], Z[M], dZ[M] = iterate(C[M], I[M], E[M], Z[M], dZ[M])
|
||||
|
||||
return E, I, Z, dZ
|
||||
return I, E, Z, dZ
|
||||
|
||||
for j in nb.prange(C.shape[1]):
|
||||
E[:,j], I[:,j], Z[:,j], dZ[:,j] = iteration(C[:,j])
|
||||
for j in numba.prange(d+1):
|
||||
I[:, j], E[:, j], Z[:, j], dZ[:, j] = iteration(S, C[:, j])
|
||||
|
||||
return E, I, Z, dZ
|
||||
return I, E, Z, dZ
|
||||
|
||||
E, I, Z, dZ = calculation(C)
|
||||
I, E, Z, dZ = iteration_numba(S, C)
|
||||
D = np.zeros(C.shape, dtype=np.float64)
|
||||
|
||||
N = abs(Z) > 2 # exterior distance estimation
|
||||
|
|
|
|||
|
|
@ -1,8 +1,85 @@
|
|||
print(
|
||||
'\n'.join(
|
||||
''.join(
|
||||
' *'[(z:=0, c:=x/50+y/50j, [z:=z*z+c for _ in range(99)], abs(z))[-1]<2]
|
||||
for x in range(-100,25)
|
||||
)
|
||||
for y in range(-50,50)
|
||||
))
|
||||
import jax
|
||||
jax.config.update("jax_enable_x64", True) # faster on GPU P100 than on GPU T4
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
import cupy as cp
|
||||
import jax.numpy as jnp
|
||||
|
||||
import decimal as dc # decimal floating point arithmetic with arbitrary precision
|
||||
dc.getcontext().prec = 80 # set precision to 80 digits (about 256 bits)
|
||||
|
||||
d, h = 100, 2000 # pixel density (= image width) and image height
|
||||
n, r = 100000, 100000.0 # number of iterations and escape radius (r > 2)
|
||||
|
||||
a = dc.Decimal("-1.256827152259138864846434197797294538253477389787308085590211144291")
|
||||
b = dc.Decimal(".37933802890364143684096784819544060002129071484943239316486643285025")
|
||||
|
||||
S = np.zeros(n+1, dtype=np.complex128)
|
||||
u, v = dc.Decimal(0), dc.Decimal(0)
|
||||
|
||||
for i in range(n+1):
|
||||
S[i] = float(u) + float(v) * 1j
|
||||
if u * u + v * v < r * r:
|
||||
u, v = u * u - v * v + a, 2 * u * v + b
|
||||
else:
|
||||
print("The reference sequence diverges within %s iterations." % i)
|
||||
break
|
||||
|
||||
x = np.linspace(0, 2, num=d+1, dtype=np.float64)
|
||||
y = np.linspace(0, 2 * h / d, num=h+1, dtype=np.float64)
|
||||
|
||||
A, B = np.meshgrid(x * np.pi, y * np.pi)
|
||||
C = (- 8.0) * np.exp((A + B * 1j) * 1j)
|
||||
|
||||
def iteration_cupy(S, C):
|
||||
|
||||
def iteration(S, C):
|
||||
I = cp.zeros(C.shape, dtype=np.intp)
|
||||
E, Z, dZ = cp.zeros_like(C), cp.zeros_like(C), cp.zeros_like(C)
|
||||
|
||||
for i in range(n):
|
||||
M = cp.absolute(Z) < cp.absolute(E) # rebase when z is closer to zero
|
||||
I, E = cp.where(M, 0, I), cp.where(M, Z, E) # reset reference orbit
|
||||
M = cp.absolute(Z) < r
|
||||
I, E = cp.where(M, I + 1, I), cp.where(M, (2 * S[I] + E) * E + C, E)
|
||||
Z, dZ = cp.where(M, S[I] + E, Z), cp.where(M, 2 * Z * dZ + 1, dZ)
|
||||
|
||||
return I, E, Z, dZ
|
||||
|
||||
I, E, Z, dZ = iteration(cp.asarray(S), cp.asarray(C))
|
||||
return I.get(), E.get(), Z.get(), dZ.get()
|
||||
|
||||
def iteration_jax(S, C):
|
||||
|
||||
def iteration(S, C):
|
||||
I = jnp.zeros(C.shape, dtype=np.intp)
|
||||
E, Z, dZ = jnp.zeros_like(C), jnp.zeros_like(C), jnp.zeros_like(C)
|
||||
|
||||
def abs2(z):
|
||||
return z.real * z.real + z.imag * z.imag
|
||||
|
||||
def iterate(i, V):
|
||||
I, E, Z, dZ = V
|
||||
M = abs2(Z) < abs2(E) # rebase when z is closer to zero
|
||||
I, E = jnp.where(M, 0, I), jnp.where(M, Z, E) # reset reference orbit
|
||||
M = abs2(Z) < abs2(r)
|
||||
I, E = jnp.where(M, I + 1, I), jnp.where(M, (2 * S[I] + E) * E + C, E)
|
||||
Z, dZ = jnp.where(M, S[I] + E, Z), jnp.where(M, 2 * Z * dZ + 1, dZ)
|
||||
return I, E, Z, dZ
|
||||
|
||||
I, E, Z, dZ = jax.lax.fori_loop(0, n, iterate, (I, E, Z, dZ), unroll=10)
|
||||
return I, E, Z, dZ
|
||||
|
||||
I, E, Z, dZ = iteration(jnp.asarray(S), jnp.asarray(C))
|
||||
return np.asarray(I), np.asarray(E), np.asarray(Z), np.asarray(dZ)
|
||||
|
||||
I, E, Z, dZ = iteration_jax(S, C) # use iteration_cupy or iteration_jax
|
||||
D = np.zeros(C.shape, dtype=np.float64)
|
||||
|
||||
N = abs(Z) > 2 # exterior distance estimation
|
||||
D[N] = np.log(abs(Z[N])) * abs(Z[N]) / abs(dZ[N])
|
||||
|
||||
plt.imshow(D.T ** 0.015, cmap=plt.cm.gist_ncar, origin="lower")
|
||||
plt.savefig("Mercator_Mandelbrot_deep_map.png", dpi=200)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue