Как нарисовать 40 соприкасающихся окружностей?

REM
REM  Программа рисует соприкасающиеся окружности
REM
REM  Формула 1 из теоремы Декарта (геометрия)
REM  k4 = k1 + k2 + k3 + 2 * SQR(k1 * k2 + k2 * k3 + k3 * k1)
REM
REM  Формула 2 из комлексной теоремы Декарта
REM  z4 = (k1 * z1 + k2 * z2 + k3 * z3 + 2 * SQR(k1*k2*z1*z2 + k2*k3*z2*z3 + k3*k1*z3*z1))/k4
REM 
REM  k1, k2, k3, k4 - вещественные числа, z1, z2, z3, z4 - комплексные числа
REM 
 
DECLARE SUB Krumm ()
DECLARE SUB Komplex ()
DECLARE SUB Index (i!, j!)
DECLARE SUB Qwrz (c0!, c1!, c00!, c11!)
 
CLS
 
CONST pi = 3.141593
DIM SHARED k(1 TO 43)
   k(1) = 1: k(2) = 1: k(3) = 1
 
DIM SHARED z(1 TO 43, 1)
   z(1, 0) = 0: z(1, 1) = 0
   z(2, 0) = 2: z(2, 1) = 0
   z(3, 0) = 1: z(3, 1) = SQR(3)
 
DIM SHARED M(4 TO 43, 3)
   M(4, 1) = 1: M(4, 2) = 2: M(4, 3) = 3
 
FOR i = 4 TO 16
   CALL Index(i, 3 * i - 8)
NEXT
 
CALL Krumm
CALL Komplex
 
SCREEN 11
WINDOW (.333, 1)-(1.667, 0)
CIRCLE (0, 0), 1, , 0, pi / 3
CIRCLE (2, 0), 1, , 2 / 3 * pi, pi
CIRCLE (1, SQR(3)), 1, , 4 / 3 * pi, 5 / 3 * pi
 
FOR i = 4 TO 43
   CIRCLE (z(i, 0), z(i, 1)), 1 / k(i)
NEXT
END
 
                ' Вычисляет номера троек соприкасающихся окружностей
SUB Index (i, j)
   M(j + 1, 1) = i
   M(j + 1, 2) = M(i, 2)
   M(j + 1, 3) = M(i, 3)
 
   M(j + 2, 1) = M(i, 1)
   M(j + 2, 2) = i
   M(j + 2, 3) = M(i, 3)
 
   M(j + 3, 1) = M(i, 1)
   M(j + 3, 2) = M(i, 2)
   M(j + 3, 3) = i
END SUB
 
                ' Вычисление декартовых координат окружностей
SUB Komplex
   FOR i = 4 TO 43
      k1 = k(M(i, 1))
      k2 = k(M(i, 2))
      k3 = k(M(i, 3))
      k4 = k(i)
      
      z10 = z(M(i, 1), 0)
      z11 = z(M(i, 1), 1)
      z20 = z(M(i, 2), 0)
      z21 = z(M(i, 2), 1)
      z30 = z(M(i, 3), 0)
      z31 = z(M(i, 3), 1)
      
      c0 = z10 * z20 - z11 * z21
      c1 = z10 * z21 + z11 * z20
      c0 = c0 * k1 * k2
      c1 = c1 * k1 * k2
 
      d0 = z20 * z30 - z21 * z31
      d1 = z20 * z31 + z21 * z30
      c0 = c0 + d0 * k2 * k3
      c1 = c1 + d1 * k2 * k3
 
      d0 = z30 * z10 - z31 * z11
      d1 = z30 * z11 + z31 * z10
      c0 = c0 + d0 * k3 * k1
      c1 = c1 + d1 * k3 * k1
 
      CALL Qwrz(c0, c1, c00, c11)
 
      z(i, 0) = (k1 * z10 + k2 * z20 + k3 * z30 + 2 * c00) / k4
      z(i, 1) = (k1 * z11 + k2 * z21 + k3 * z31 + 2 * c11) / k4
   NEXT
END SUB
 
                ' Вычисление кривизны окружностей
SUB Krumm
   FOR i = 4 TO 43
      k1 = k(M(i, 1))
      k2 = k(M(i, 2))
      k3 = k(M(i, 3))
      k(i) = k1 + k2 + k3 + 2 * SQR(k1 * k2 + k2 * k3 + k3 * k1)
   NEXT
END SUB
 
              ' Извлечение квадратного корня из числа a + bi
SUB Qwrz (a!, b!, c00!, c11!)
    c11 = SQR((SQR(a ^ 2 + b ^ 2) - a) / 2)
    c00 = b / (2 * c11)
END SUB