from math import sin, cos
def shrink_rotate(x, a, coef=1., angle=0.):
    c = x.real
    d = x.imag
    e = a.real
    f = a.imag
    axx = (x.real - a.real) * coef
    axy = (x.imag - a.imag) * coef
    x_new = axx * cos(angle) - axy * sin(angle)
    y_new = axx * sin(angle) + axy * cos(angle)

    return x_new + y_new  * 1j

def shrink_rotate_conj(x, a, coef=1., angle=0.):
    c = x.real
    d = x.imag
    e = a.real
    f = a.imag
    axx = (x.real - a.real) * coef
    axy = (x.imag - a.imag) * coef
    x_new = axx * cos(angle) - axy * sin(angle)
    y_new = axx * sin(angle) + axy * cos(angle)

    if y_new < f:
        y_new += 2 * abs(f - y_new)
    else:
        y_new -= 2 * abs(f - y_new)

    return x_new + y_new  * 1j

z = 0.5 + 0.*1j
max_iter = 100000
funcs = [
    (lambda t: shrink_rotate(t, 0. + 1.*1j, coef=0.5, angle=0.)),
    (lambda t: shrink_rotate(t, 1. + 0.*1j, coef=0.5, angle=0.)),
    (lambda t: shrink_rotate(t, -1. + 0.*1j, coef=0.5, angle=0.))
]

res = np.zeros((1001, 1001), dtype=int)

for n_iter in range(max_iter):
    n_func = np.random.choice(len(funcs))
    z = funcs[n_func](z)
    x = int(z.real * 500 + 500)
    y = int(z.imag * 500 + 500)
    res[1000 - y][x] = 1

plt.figure(figsize=(20, 20))
plt.imshow(res, cmap='gray')