https://twitter.com/MauriceMeilleur/status/1338222143503527943

noMagic.py

def lerp(value1, value2, factor):
    return value1 + (value2 - value1) * factor

size(400, 400)

R = (sqrt(5) - 1) / 2


def F(x1, y1, r, s, t, n):
    x2 = x1 + r * cos(s)
    y2 = y1 + r * sin(s)    
    x3 = x1 + r * cos(t)
    y3 = y1 + r * sin(t)  
    q = t - s

    if r < 50 * R ** n:
        polygon((x1, -y1), (x2, -y2), (x3, -y3), close=False)
    else:
        if n:
            F(x3, y3, r, t + q * 2, t + q * 5, 0)
            F(lerp(x3, x2, R), lerp(y3, y2, R), r * R, t + q, t + q * 2, 1)
        else:
            F(x2, y2, r * R, t + q * 2, t + q * 3, 0)
            F(x2, y2, r * R, t + q * 4, t + q * 3, 0)
            F(lerp(x2, x1, R), lerp(y2, y1, R), r * (1 - R), t + q, t + q * 4, 1)

stroke(0)
fill(None)
translate(0, height())
F(-800, 0, 1600, pi/5, 0, 0)