跳转到内容

File:Exponential Collatz Fractal.jpg

页面内容不支持其他语言。
這個文件來自維基共享資源
维基百科,自由的百科全书

原始文件 (30,720 × 17,280像素,文件大小:41.13 MB,MIME类型:image/jpeg


摘要

警告 部分浏览器在浏览此图片的完整大小时可能会遇到困难:该图片中有数量巨大的像素点,可能无法完全载入或者导致您的浏览器停止响应。 交互式大图查看器
描述
English: A Collatz fractal for the interpolating function . The center of the image is and the real part goes from to .
日期
来源 自己的作品
作者 Hugo Spinelli
其他版本
new file 本图像是 File: Exponential Collatz Fractal.png 中原始图像 PNG 的 JPEG 版本。

一般在显示维基共享资源中的文件时应使用本 JPEG 版本,以减小需传输图片的大小。 不过,编辑时应使用原始 PNG 版本以避免叠代丢失,且需同时更新两种版本。请勿编辑本版本 请参阅此处了解更多信息。

Source code
InfoField
Python code
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

许可协议

我,本作品著作权人,特此采用以下许可协议发表本作品:
Creative Commons CC-Zero 本作品采用知识共享CC0 1.0 通用公有领域贡献许可协议授权。
采用本宣告发表本作品的人,已在法律允许的范围内,通过在全世界放弃其对本作品拥有的著作权法规定的所有权利(包括所有相关权利),将本作品贡献至公有领域。您可以复制、修改、传播和表演本作品,将其用于商业目的,无需要求授权。

说明

添加一行文字以描述该文件所表现的内容
An exponential Collatz fractal with smooth coloring based on divergence speed.

此文件中描述的项目

描繪內容

文件历史

点击某个日期/时间查看对应时刻的文件。

日期/时间缩⁠略⁠图大小用户备注
当前2023年10月6日 (五) 19:372023年10月6日 (五) 19:37版本的缩略图30,720 × 17,280(41.13 MB)Hugo SpinelliUploaded own work with UploadWizard

以下页面使用本文件:

全域文件用途

以下其他wiki使用此文件: