сохраню ссылку на семинар по экспериментальной математике, который Zeilberger с коллегами проводит уже больше 20 лет

https://sites.math.rutgers.edu/~zeilberg/expmath/

есть и видеозаписи

суммы двух квадратов обсуждают часто, а родственные утверждения — реже

для этого есть объективные причины, но все же попробуем еще что-то… вот, например, таблица чисел вида X²+5Y²

как и в прошлый раз, посмотрим на представимость нечетных чисел — и сразу возникают какие-то идеи и какие-то вопросы, правда?

задача для начала знакомства с темой: доказать, что если числа M и N представимы, то и число MN представимо

более искушенным читателям предлагаю обратить внимание на то, что числа 3, 23, 43 не представимы в нужном виде, хотя кратные им числа нетривиальным образом представимы (для X²+Y² так не бывает)

***

в экселе ставить простые эксперименты легко, и выглядит сразу нагляднее, чем код

но все ж для данного d поискать в представимости N в виде X²+dY² период, да еще не учитывая при этом составные N, удобнее на питоне (код будет в комментариях)

***

Эйлер размышлял над этим сюжетом примерно 40 лет и открыл много всего интресеного, включая квадратичный закон взаимности

поразительно, что Эйлер нашел (хотя и не доказал) правильный ответ и, например, для d=27, где… ну можете сами запустить программу и посмотреть

вокруг этого сюжета есть замечательная книжка Cox'а «Primes of the form x²+ny²» (она быстро становится не очень простой — по крайней мере, я теряюсь — но первая глава заинтересованным старшеклассникам вполне доступна)

как делить числа?

пусть у нас есть числа a и b и мы хотим быстро посчитать a/b с большой точностью (а складывать-умножать числа мы уже умеем)

можно сдвинуть числитель и знаменатель на степень двойки, так что достаточно научиться находить 1/b для b между, скажем, 1/2 и 1

в этом посте нет экспериментальной математики, но так понравилась история, что не могу не поделиться — спасибо рассказавшему про это А.Гасникову (все ошибки, естественно, на моей совести)

1.

когда надо вычислить значение функции, у меня первый рефлекс — разложить ее в ряд Тейлора, т.е. в данном случае просто воспользоваться бесконечной геометрической прогрессией:

если b=1-q, то 1/b=1+q+q²+q³+… — сиди и вычисляй столько членов, сколько тебе нужно (так как |q|<1/2, рано или поздно всё получится)… но сколько нужно? чтобы найти N (двоичных) знаков после запятой нужно взять ~N членов, т.е. сделать ~N умножений, и это не очень вдохновляет

в конце концов, уравнение f(x)=0 для монотонной функции b всегда можно решать методом деления пополам (выбираем ту половину, на концах которой у f разные знаки, повторяем процесс) — уже это позволяет искать N знаков после запятой за ~N действий (для деления такой способ так же известен как деление в столбик)

но оказывается, что делить можно и намного быстрее, чем в столбик!

2.

если функция достаточно хорошая, то уравнение f(x)=0 можно быстро приближенно решать при помощи метода Ньютона

(
напомню идею: если x — приближенное значение корня, то рядом ним график функции недалеко ушел от касательной, поэтому в качестве следующего приближения можно взять пересечения касательной с нулем, т.е. x→x-f(x)/f'(x)

и в некоторой окрестности корня метод Ньютона, если производная в этом корне не равна 0, сходится очень быстро: за итерацию погрешность ~возводится в квадрат
)

на первый взгляд, нам метод Ньютона не поможет, так как в него входит деление — но тут происходит чудо: если сформулировать нашу задачу как задачу поиска нуля функции f(x)=1/x-b, то в методе Ньютона все деления сокращаются: f/f'=(1/x-b)⋅(-x²)=-x⋅(1-bx)

и получается рецепт x→x⋅(2-bx), который позволяет получить N знаков числа 1/b всего за ~log(N) операций (за каждую операцию количество верных знаков ~удваивается)

можно проверить, как это работает:

from mpmath import *
mp.dps = 300

b = mpf(57)/100
x = mpf(1)
print("1/"+nstr(b,30))
print(0,":",nstr(x,80))
for k in range(10):
    x = x*(2-b*x)
    print(k+1,":",nstr(x,80),
          "diff:",nstr(abs(1/b-x),2,min_fixed=1))
print("T :",nstr(1/b,80))

3.

что это всё-таки за странная формула x→x(2-bx), можно ли это связать с чем-то более знакомым?

оказывается, это просто способ быстро вычислять частичные суммы всё той же геометрической прогрессии! действительно, если b = 1-q, и x = (1-q^N)/(1-q), то x’ = (1-q^N)/(1-q) ⋅ (1+q^N) = (1-q^{2N})/(1-q) — т.е. за шаг мы переходим от суммы первых N членов к сумме первых 2N членов геометрической прогрессии — немножко похоже на быстрое возведение в степень, только еще формулы сокращенного умножения знать надо

===

что можно ускорять какие-то базовые операции, меня впечатляет; вот небольшой текст про это в Мат. составляющей (в т.ч. про быстрое умножение): https://book.etudes.ru/articles/fast-arithmetic/

метод Ньютона здесь уже появлялся раньше, и будет, думаю, обсуждаться еще

дайджест комментариев: разбиения на доминошки

коллеги, спасибо большое за содержательные комментарии и вообще

1.
С.Шашков сделал веб-версию перемешивания ацтекского брильянта: https://shashkovs.ru/i/Aztec.html

Нравится и как конкретно это выглядит, и вообще идея превращения таких программ в веб-страницы — особ. удобно, если хочется поделиться на кружке, докладе и т.п.

Переход от несложного питона к джаваскрипту выглядит посильным — мб попробую при случае что-то сделать и на джаваскрипте.

2.
Л.Петров обращает внимание на то, что случайное разбиение ацтекского брильянта на доминошки радикально быстрее генерировать не применяя много случайных флипов, а при помощи domino shuffling.

В качестве популярного введения — советую видео https://youtu.be/Yy7Q8IWNfHM Mathologer'а. Для тех, кто читал брошюру Е.Смирнова про ацтекские брильянты, — это примерно то же, что описанное там «расширение площадей».

3.
Р.Гусарев напоминает, что разбиения квадрата 8×8 на доминошки намного быстрее считать не в лоб, а «динамикой по профилю».

Эту идею давайте в таком виде упакую. Если считать просто разбиения прямоугольника, say, 3×N, то эти числа P(n) никакой очевидной рекурренте не удовлетворяют. Почему? Ну просто потому, что если мы выкидываем все доминошки, покрывающие последний столбец, то остается не прямоугольник, а прямоугольник с дырками в самом правом столбце. Но это значит, что если думать про тройку (P(n),Q(n),P(n-1),Q(n-1)), где Q(n) — количество разбиений прямоугольника 3×N без верхней правой клетки, P(n-1) — количество разбиений без всего правого столбца, S(n-1) — только с одной клеткой в самом правом столбце, то следующая четверка линейно выражается через предыдущую!

Реализовал это вот так (и теперь, действительно, даже разбиения прямоугольника 8×64 считаются мгновенно):

def is_good(mask,n):
    # mask кодирует последовательность из n нулей и единиц
    # функция проверяет, можно ли замостить нули доминошками
    if mask == 0:
        return n%2 == 0
    if mask%4 == 0:
        return is_good(mask>>2,n-2)
    if mask%4 == 2:
        return False
    return is_good(mask>>1,n-1)

def tilings(n,m):
    ext = [ [1 if mask&perp==0 and is_good(mask+perp,n) else 0
            for perp in range(2**n)] for mask in range(2**n)]
    # ext[mask][perp]: можно ли положить перпендикулярные
    #   нашему ряду доминошки, чтобы не задеть маску,
    #   а остаток чтобы разбился на доминошки в ряду
    ans = [1 if mask==0 else 0 for mask in range(2**n)]
    for _ in range(m):
        newans = [0] * (2**n)
        for mask in range(2**n):
            for oldmask in range(2**n):
                newans[mask] += ext[mask][oldmask]*ans[oldmask]
                # видно, что шаг есть умножение матрицы на вектор
        ans = newans
    return ans[0]

print(tilings(8,64))

Дима Швецов напомнил про такой сюжет

Квантик пытается экспериментально изучить сумму обратных простых

Ноутик воспользовался освоенной им недавно библиотекой matplotlib и построил график частичных сумм до разных N (см. выше), нашел что у правого конца графика принимаются значения между 2.88732 и 2.88733, и теперь пытается составить из e и π выражение с таким примерно значением

а чему на самом деле равна эта бесконечная сумма?

...

ряд расходится — просто медленно (эмпирическое объяснение: случайное число ~n простое с вероятностью ~1/log(n), так что сумма растет примерно как интеграл от 1/(x log(x)), т.е. примерно как log log(N))

так что компьютерные эксперименты — это классно, но мало помогают, если нет общего математичекого понимания, что происходит

it is the theory which decides what can be observed

что делать с ответами?

начнем давайте с модельной задачи… вот хотя бы с теми же доминошками: посчитаем количество P(N) разбиений на доминошки прямоугольника 3×N

как объяснялось выше, на эти числа [вместе с количествами разбиений прямоугольников без угла] можно написать линейную рекурренту (она получается из того, что у правого края либо все доминошки горизонтальные, либо одна вертикальная), поэтому программирование не требуется и можно считать всё в экселе

хорошо, сделали файл, в котором посчитаны сколько-то первых значений P(N) — и что дальше?

можно построить график (см. картинку над постом), но на нем мало что видно: начало графика кажется совершенно горизонтальным, в конце резкий рост

выбранное представление данных пока не позволяет легко увидеть, что именно происходит

когда целые числа быстро растут, то совсем грубо можно смотреть просто на количество цифр в P(N) — ну или если не на таком детском языке, то можно говорить про логарифмическую шкалу для оси y

если это сделать и прострить график не P(2N), а log P(2N), то… ну это практически как маленькое чудо выглядит (и на школьников, вроде, производит некоторое впечатление)

продолжение следует

что делать с ответами? (продолжение)

если построить график не количеств разбиений прямоугольников 3×2N на доминошки, а логарифмов этих количеств, то асимптотика сразу видна: график визуально неотличим от прямой

что это конкретно за прямая? можно не пытаться подобрать коэффициенты на глаз и т.п., а воспользоваться функцией экселя LINEST

если log P(2N) ≈ kN+b, то P(2N) ≈ c⋅λ^N (для λ=10^N, c=10^b)

хорошее ли получается приближение? вот фиолетовым цветом внизу таблицы показаны эти самые приближения c⋅λ^N, а рыжим цветом в последней строчке — величина ошибки

видно, что приближение суперхорошее, дальше можно прямо угадать точную формулу для P(2N) (указание: на что похоже приближенно найденное λ? — надо бы, кстати, написать программу, которая на подобные вопросы отвечает)

эксель в комментариях

===

конечно, настолько хорошо всё видно не всегда — дальше планируется не настолько модельный пример

различаются ли вычисленные приближенно рациональные и иррациональные числа?

про математику здесь хорошо написано в самом начале книги «Математический дивертисмент» Табачникова и Фукса (которую вообще всячески рекомендую)

а я только покажу программу, которая различает рациональные и иррациональные числа: по числу, например, 1,414213562 говорит, что оно иррациональое, а по числу 1,470588235 сразу понимает, что это 25/17

import math
from fractions import Fraction

def contfrac(x,n):
    ans = [0]*n
    for i in range(n):
        ans[i] = math.floor(x)
        if x-ans[i] == 0:
            return ans
        else:
            x = 1/(x-ans[i])
    return ans

def contvalue(arr):
    if len(arr)==1 or arr[1]==0:
        return Fraction(arr[0])
    return arr[0]+1/contvalue(arr[1:])

x = 1.470588235

digits = contfrac(x,17)
ans = None
for i in range(1,len(digits)):
    if digits[i]>=1000 or digits[i]==0:
        ans = contvalue(digits[:i])
        break
if ans is None:
    print("irrational")
else:
    print(ans)

хотел еще сделать, чтобы распознавались квадратичные иррациональности (и можно было ввести λ из прошлого поста, например), но пока руки не дошли

на картинке сверху — тождества¹ из заметки С.Маркелова в Мат. просвещении, и там предлагается придумать обобщения

¹ там только есть опечатка… найдите

программа в комментариях — говорит, суммы каких косинусов надо взять для произвольного p вида 3k+1, а также какому кубическому уравнению они удовлетворяют (и на всякий случай численно проверяет, удовлетворяют ли)

(upd) а также находит формулу для суммы S кубических корней из этих сумму косинусов, шоб было совсем как в заметке

p: 13
primitive root: 2
partition of cosines: [3, 11] [7, 9] [1, 5]
values of trigsums: -0.136945 -0.688601 1.325547
cubic polynomial: 8t³-4t²-8t-1
P(trigsums): -0.0 0.0 0.0
S³ = (3³√-13+7)/2

p: 73
primitive root: 5
partition of cosines: [13, 19, 25, 29, 31, 39, 53, 55, 57, 59, 67, 71] [1, 3, 7, 9
, 17, 21, 27, 43, 49, 51, 63, 65] [5, 11, 15, 23, 33, 35, 37, 41, 45, 47, 61, 69]
values of trigsums: -2.475085 2.40906 0.566026
cubic polynomial: 8t³-4t²-48t+27
P(trigsums): 0.0 0.0 -0.0
S³ = (3³√219-17)/2

теорема Рамануджана о том, как посчитать сумму кубических корней из корней данного кубического уравнения, обсуждается например в Кванте

суммы косинусов и кубические уравнения

сегодня кода не будет, а будет математический комментарий к предыдущему посту (и не совсем популярный, сорри)

1.
напомню для начала квадратичный случай

пусть P простое число… и пускай вида 4k+1

квадратичная сумма Гаусса — это способ выразить √P через корни степени P из единицы с разными знаками (если показатель степени — квадратичный вычет mod P, то берем знак +1, для невычета — -1)

можно записать его в виде ∑cos(n² 2π/P)=√P

(например, для p=5: cos(2π/5)=(√5-1)/4, а вся сумма состоит из 4 равных косинусов и единицы)

2.
пускай теперь P=3k+1, а нас интересует тригонометрическая сумма G=∑cos(n³ 2π/P)

Галуа говорит нам, что G — корень кубического уравнения с рациональными коэффициентами, но какого конкретно?

полезно начать не с суммы G, а с кубической суммы Гаусса g

определяется она аналогично квадратичному случаю, но знак +1 теперь получают степени, являющиеся кубическими вычетами, еще два класса — знаки ω и ω² (ω — комплексный кубический корень из единицы)

(другими словами, g — линейная комбинация трех корней искомого уравнения, на которой группа Галуа действует просто умножением на ω)

снова |g|=√P, но теперь что-то хорошее получается при возведении не в квадрат, а в куб: g³=PП, где ПП' — разложение p на простые в Z[ω], а штрих обозначает комплексное сопряжение

3.
вернемся к тригонометрической сумме

G=g+g', поэтому (пользуясь описанными выше свойствами g)

G³=g³+g'³+3gg'(g+g')=P(П+П')+3PG,

т.е. G корень кубического уравнения x³-3Px-AP, где A=2ReП

чтобы получить уравнение буквально на суммы, которые рассматривались в предыдущем посте, нужно еще сделать замену x=-3t+1

4.
если хочется вычислить A без использования арифметики в Z[ω], то можно воспользоваться тем, что если P=ПП' в Z[ω], A=2ReП, то

4P=A²+27B²

такое представление единственно… с точностью до знаков; “наше” A — это то, которое дает остаткок 1 mod 3 (это вопрос про знак гауссовой суммы, который выше был заметен под ковер… там пусть и остается)

можно описать A и по-другому: количество решений по модулю P уравнения X³+Y³=1 равно P-2+A

5.
дискриминант уравнения x³-3Px-AP=0 как раз равен (27PB)²

вообще через корни из единицы (или, если угодно, через триг. функции рац. аргумента) можно выразить решения тех и только тех неприводимых кубических уравнений с рациональными коэффициентами, у которых дискриминант является полным квадратом (аналог для произвольной степени: группа Галуа должна быть абелевой)

уравнение x³=P, например, не подходит, а вот наше — вполне себе

надеюсь, что в написанном выше можно уже увидеть контуры рецепта для решения «в косинусах» произвольного кубического уравнения с квадратным дискриминантом

6.
если P=9m²-3m+1, то 4P=(3m-2)²+27m². В этом случае свободный член уравнения на триг. суммы является полным кубом (числа m, собственно) и формула Рамануджана для суммы кубических корней из корней уравнения дает особенно простой ответ:

2S³ = 3³√(mP)-6m+1

например, для P=73 имеем m=3 и ответ ³√((3³√219-17)/2); для P=31 имеем m=-2 и опечатку в заметке в Мат. просвещении

L-системы, или Как нарисовать снежинку и дерево

1.
начнем со слова A и дальше будем итерировать такие замены: A→AB, B→A (на каждой итерации одновременно заменяются все буквы по указанным правилам)

будут получаться всё большие отрезки интересного слова Фибоначчи — про него можно прочитать в статье Вити Клепцына в «Квантике» (№№5-6 за 2020 год)

такого рода правила (“L-системы”) легко реализовать — например, так:

def lsystem(seed,rules,iterations):
    for _ in range(iterations):
        seed = "".join(
            rules[letter] if letter in rules else letter
            for letter in seed)
    return seed

print(lsystem("A",{"A":"AB","B":"A"},5))

2.
такого рода правила порождают слова с фрактальной структурой… так что можно попробовать превратить их во фрактальные картинки

можно начать, например, со снежинки Коха: мы стартуем с правильного треугольника, а дальше на каждой итерации заменяем каждый отрезок на зубец («_ → _/\_»)

это можно закодировать так:

snowflake = lsystem("F++F++F",{"F":"F-F++F-F"},4)

где “F” означает «идем вперед и рисуем отрезок», а “+” и “-” — изменения направления

(про такое, кстати, тоже писали в «Квантике»: Валя Кириченко и Владлен Тиморин объясняли, как это связано со словом Туэ-Морса и проч.)

осталось это реально нарисовать — например, при помощи matplotlib:

import math
import matplotlib.pyplot as plt

def turtleprint(string,alpha0=0,alpha=90,step=1):
    fig = plt.figure()
    fig.set_size_inches(6,6)
    ax = fig.gca()
    ax.set_aspect('equal')

    x,y = 0,0
    heading = alpha0
    for command in string:
        if command == "F":
            x0,y0 = x,y
            x = x0 + step*math.cos(math.radians(heading))
            y = y0 + step*math.sin(math.radians(heading))
            plt.plot([x0,x],[y0,y],color='k')
        if command == "+":
            heading += alpha
        if command == "-":
            heading -= alpha

    plt.axis('off')
    plt.show()

для 4 итераций получается буквально картинка сверху

можно в таком духе генерировать и, скажем, несамопересекающиеся кривые, приближающиеся к заполнению квадрата — например, вот так:

turtleprint(lsystem("-L",{"L":"LF+RFR+FL-F-LFLFL-FRFR+","R":"-LFLF+RFRFR+F+RF-LFL-FR"},4))

(продолжение следует)

L-системы, или Как нарисовать снежинку и дерево

(продолжение, начало выше)

биолог Аристид Линденмайер в конце 1960-х придумал L-системы именно как модель для роста растений — и невозможно не написать про это два слова

сразу скажу, что есть отличная книжка Prusinkiewicz, Lindenmayer. The Algorithmic Beauty of Plants, http://algorithmicbotany.org/papers/#abop (спасибо Никите С. и Роме Г., советовавшим эту книжку) — и там все подробно step by step объясняется, лучшее мб вместо поста ее и читать )

===

3.
в предыдущем посте мы генерировали строки и смотрели на них как на описания путей (символ F означает «двигаться вперед», повороты кодируются символами + и -)

если хочется рисовать деревья, то полезно добавить какую-то возможность ветвления

для этого добавим такое правило для обработки символов [ и ]: после рисования фрагмента в квадратных скобках «черепашка» возвращается туда, где была в его начале

технически для этого в функцию turtleprint добавим

stack = []
…
if command == "[":
    stack.append((x,y,heading))
if command == "]":
    x,y,heading = stack.pop()

тогда дерево можно нарисовать в духе

seed="F"
rules={"F": "F[+F][-F]F"}

(стартуем с отрезка, на каждом шаге заменяем отрезок на цепочку из двух с растущими из середины в обе стороны ветками; на рисунке показано всего 2 шага)

получается, правда, такое совсем уж математическое дерево… чтобы сделать его более растительным, можно сделать два вида “клеток” — "F" пусть будет просто расти, а из "X" пусть происходит еще ветвление в разные стороны

seed="X"
rules={"F": "FF", "X": "F[+X][-X]FX"}

уже не так плохо, но слишком симметрично… перекосить можно так, скажем:

seed="X"
rules={"F": "FF", "X": "F-[[X]+X]+F[+FX]-X"}

(код целиком положу в комментарии)

***

(не слишком сложными) L-системами можно пользоваться и прямо в техе (в tikz есть для этого бибиотека):

\usetikzlibrary{lindenmayersystems}
\begin{tikzpicture}[rotate=65]
\draw[green!60!black] l-system[l-system={axiom=F,rule set={F -> F[+F]F[-F]}, order=4, angle=25,step=3pt}];
\end{tikzpicture}

но в питоне, конечно, легче менять любые правила игры по своему желанию

===

вижу иногда, как люди постят generative art, где код умещается примерно в один твит, а на картинках облака летят, деревья растут, вселенные встают и рушатся, и всё такое

здесь вещи довольно примитивные, но всё-таки от того, что вот такая одна строчка фактически рисует узнаваемое дерево (ну хорошо, куст), и вообще от темы в целом — такого же рода вайбы

Посмотрим на то, для каких простых P уравнение X²=D разрешимо в остатках по модулю P.

На компьютерных экспериментах хорошо видно, что возникает удивительная периодичность — вот, например, для D=7 (для непростых модулей вместо ответа печатается точка):

x²=7 mod p
apparent (weak) period is 28
..++.-.+...-.-...-.+...-....
.+.+.....+...-.-...+.....+..
...+.-.....-...-.-.....-...+
.....-.......-...-.+...-.+..
.+.............-...+.....+.+
.........+.-.....-.....-...+
.....-.....-.-.........-.+..
.+.+...........-...........+
...+.-...+.....-.-.........+
.....-.....-.....-.+.....+..

Ответ всегда зависит только от остатка P mod 4D — это квадратичный закон взаимности (в форме, открытой Эйлером).

Можно было бы надеяться, что похожая история будет и в кубическом случае… но нет, ничего похожего. Простые P, для которых уравнение X³=2 (например) разрешимо mod P, выглядят хаотично, никакого периода нет. Единственная относительно видная закономерность — что доля простых, для которых X²=D разрешимо, всегда (ну если D не куб ) примерно 2/3.

x³=2 mod p
no apparent period
..++.+.-...+.-...+.-...+.....+.+.....-...+.+...+.....+.....+.-.....-...+.-.....-...+.....+.......-...+.-...+.+...+.............+...+.....+.-.........+.-.....+.....-...+.....+.....+.-.........+.-...+.-...........-...........+...+.+...+.....+.-.........+.....+.....+.....+.-.....+...+.+.........+.............(…)
+: 112 -: 56 rat: 0.67

Но можно посмотреть и на другие кубические многочлены — и иногда ответ снова определяется остатком P по какому-то модулю. Например:

x³-9x-9=0 mod p
apparent (weak) period is 9
..-+.-.-.
..-.-...+
.+...-...
..-.-....
.+...-.-.
..-.....+
.....-.-.
....-...+
.+.....-.
..-.....+
.......-.
.
+: 8 -: 17 rat: 0.32

x³-x²-4x-1=0 mod p
apparent (weak) period is 13
..--.+.-...-.
+...-.-...-..
...-.+.....-.
..-.-...+....
.+.....-.-...
..-...-.+....
.+...+.....-.
......-...-.+
...-.+...-...
..........-..
.+.....-.-...
......-.+....
.+.....-...-.
....-.....-.+
.........-.-.
..-.-
+: 14 -: 32 rat: 0.3

Можно подумать, чем одни кубические многочлены отличаются от других (если гадать станет невыносимо, то посчитайте дискриминанты).

===

Это всё можно и в экселе делать — но в комментарии положу небольшую программу на питоне.

А для математического контекста — пусть здесь будет популярная лекция Фрэнка Калегари на ICM-2022: https://youtu.be/EDsK-8SBx-g (там как раз с таких вопросов все и начинается)

справа можно видеть фрагмент квазипериодического замощения плоскости

в нем участвуют равнобедренные треугольники с углами при вершине π/5 (красные, «A») и 3π/5 (синие, «B»)

они замечательны тем, что A можно разбить на уменьшенные копии A,B,A, ну а B можно разбить на уменьшенные копии A,B — и если начать с А и итерировать такие замены, то можно думать, что мы собираем из треугольников A и B всё большую копию треугольника¹ A (в левой половинке картинки — первая пара итераций)

такая мозаика — одна из вещей, про которые при создании канала думал, что хорошо бы ее нарисовать, но не очень понятно как

а вечером подумал, что это просто L-система — только параметрическая: кроме буквы A/B нужно помнить, как именно треугольник расположен на плоскости (и правила замены эти параметры должны правильно менять) — так что можно быстренько реализовать

¹ а чтобы получить замощение плоскости, можно, скажем, стартовать с 10 треугольников A с общей вершиной



положение треугольника решил хранить в виде пары комплексных чисел² (преобразования z→az+b, переводящего эталонный треугольник в наш) и написал такой шаг для получающейся параметрической l-системы:

phi = (math.sqrt(5)+1)/2
rot = math.cos(math.pi/5)+math.sin(math.pi/5)*1j

def step(state):
    for atom in state:
        c, a, b = atom
        if c=='A':
            yield ('A',a,b*phi)
            yield ('B',a*(rot**4),(a+b)*phi)
            yield ('A',a*(rot**3),(a+b)*phi)
        if c=='B':
            yield ('A',a,b*phi)
            yield ('B',a*(rot**4),(a+b)*phi)

state = [('A',rot**i,0+0j) for i in range(10)]
for _ in range(6):
    state = step(state)

по сути на этом все! — остается только дописать код для рисования треугольничков… ну программа целиком будет в комментариях

² уже засомневался, так ли это удачно — потому что для настоящей мозаики Пенроуза треугольники полезно и переворачивать