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

в нем участвуют равнобедренные треугольники с углами при вершине π/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)

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

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

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

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

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

попросил chatgpt сконвертировать реализацию RSK из прошлого поста в javascript, поправил минимально оформление — и пожалуйста, вот веб-версия:

https://dev.mccme.ru/~merzon/compmath/RSK.html

если для экспериментов нужен повод — поучительно посмотреть, скажем, что происходит при этом соответствии с инволюциями

будем считать количество трехмерных диаграмм Юнга («фигур из кубиков в углу комнаты», между кубиком и каждой из трех стенок не должно быть пустот) данного объема

их оказывается столько же, сколько матриц из целых неотрицательных чисел, если число в клетке (i, j) считать с весом i+j-1 (так что, например, производящая функция для всех трехмерных диаграмм равна Π(1-q^i)^{-i})

это может показаться смутно гармоничным (если в плоском разбиении на клетке (i,j) лежит кубик, то и на клетках (i',j) и (i,j') при всех i'<i и j'<j лежат кубики… и всего таких кубиков минимум i+j-1…), но если пытаться придумать конкретную биекцию, возникает ощущение магии

это тоже магия соответствия RSK


выше была версия RSK, переводящая перестановку в пару стандартных таблиц одинаковой формы… здесь используется версия, переводящая произвольную матрицу из целых неотрицательных чисел в пару полустандартных таблиц одинаковой формы — а если рассечь трехмерную диаграмму Юнга по диагонали то как раз получится обычная диаграмма Юнга, в клетки которой можно записать, сколько кубкиков растет в одну сторону, а можно — сколько в другую

реализация — в комментариях


дальше решил попросить chatgpt сделать веб-версию с рисованием кубиков — и это была прямо суперразочаровывающая трата времени (все делает неправильно, не способен реагировать на замечания, да еще при доработке рисования кубиков все время норовит удалить собственно RSK)

плюнул и за 10 минут написал matplotlib-рисование кубиков (результат на картинке)… ну может еще и сделаю веб-версию, не знаю

***

у нас с Мишей Берштейном есть текст в МатПросвещении, где рассказывается в том числе про подсчет трехмерных диаграмм Юнга — но там никакого RSK как раз нет, другое есть: https://www.mathnet.ru/rus/mp823

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

сколькими способами число N можно представить в виде суммы?

(разбиения 1+2 и 2+1 и т.п. считаются одинаковыми)

начать эксперименты можно и с примитивного перебора — например, такого:

def partitions(n,l=None):
    # все разбиения числа n
    # (на числа не превосходяшие l)
    if l is None: l=n
    if n==0: yield []
    if n>0:
        for i in range(l,0,-1):
            # рекурсивный перебор по максимальному элементу
            yield from map(lambda partition: [i]+partition,
                           partitions(n-i,i))

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

можно посмотреть на количества таких разбиений для разных N — на картинке Женя Смирнов с ассистентами использует умный способ это считать — а вот, наоборот, предельно прямолинейный (только всё же заменим рекурсию на динамику, чтобы работало и для относительно больших N):

def count_partitions_upto(N):
    partitions = [[0 for _ in range(N+1)] for _ in range(N+1)]
    for l in range(N+1):
        partitions[0][l] = 1
    for n in range(1,N+1):
        for l in range(1,N+1):
            for i in range(1,l+1):
                partitions[n][l] += partitions[n-i][i]
    return partitions

получаются какие-то числа… что дальше с ними делать? как обсуждалось выше, полезно бывает посмотреть на график — и, так как p(N) быстро растет, лучше смотреть на график log p(N)

картинки (и еще пара слов) в комментариях

дальше мб будет что-то про p(N) по разным модулям

в дополнение к предыдущему — количества p(n) разбиений — по разным модулям

ну… заметить кое-что можно, но не так-то просто (в частности, если кто-то надеялся на периодичность — то увы)

спойлер:
p(5n+4) делится на 5
p(7n+5) делится на 7
p(11n+6) делится на 11

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

в качестве мат. контекста можно прочитать короткую заметку Ramanujan’s congruences and Dyson’s crank (G.Andrews, K.Ono)

***

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

правда думал заодно и саму p(n) в экселе посчитать — но пару раз запутался в формулах и бросил (тут, конечно, с питоном намного проще) — и числа p(n) в таблицу скопировал готовые

***

заодно решил посмотреть на количества разбиений с разными рангами (контекст e.g. в заметке выше) — вот код для этого:

def partitions_upto(N):
    partitions = [[] for _ in range(N+1)]
    partitions[0] = [(".",0,0)]
    for n in range(1,N+1):
        for i in range(n,0,-1):
            for pstr,pmax,plen in partitions[n-i]:
                if i>=pmax:
                    partitions[n].append((str(i) if pstr=="." else f"{str(i)}+{pstr}", i, plen+1))
    return partitions

mod = 7
k = 4
partitions = partitions_upto(mod*k)
print(f"mod = {mod}")
for n in range(mod*k+1):
    ans = [0 for _ in range(mod)]
    for pstr,pmax,plen in partitions[n]:
        ans[(pmax-plen) % mod] += 1
    print(f"{n} — {len(partitions[n])} ({len(partitions[n]) % mod}):", *ans)

о том, как (не) стоит считать π, и о магии им. Эйлера


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

если взять arctg(x) = x - x³/3 + x⁵/5 - …, подставить¹ x=1, то получается замечательно простая формула Лейбница π/4 = 1 - 1/3 + 1/5 - 1/7 + … — чего еще желать?..

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

в прошлый раз для сотни правильных знаков хватило 4 итераций, сейчас попробуем взять… ну, скажем, 1000 членов:

from mpmath import *

n = 1000
mypi = mpf(0)
s = mpf(4)
for k in range(n):
    mypi += s / mpf(2*k+1)
    s = -s

print(n,nstr(mypi,10),
      "diff:",nstr(pi-mypi,3,min_fixed=1))

1000 3.140592654 diff: 1.0e-3

видно, что погрешность ~1/1000… и вообще, можно понять, что для N слагаемых такой способ дает погрешность ~1/N — то есть чтобы получить 100 знаков нам нужно сложить ~10^(100) слагаемых (!!!)


2.
кажется, что это совсем тупиковая ветвь… и все же, не будем унывать сразу

посмотрим, скажем, на сумму 5 млн членов

погрешность получается ожидаемая, порядка 1/N, т.е. шесть цифр после запятой верны, а седьмая уже нет (получается 3,1415924… вместо 3,1415926…)

но то, что видно дальше, иначе как магией не назовешь: например, следующая дюжина цифр снова правильная!

from mpmath import *
mp.dps = 100

n = 2500000
mypi = mpf(0)
for k in range(1,n+1):
    mypi += mpf(8) / ((4*k-3)*(4*k-1))

for (s1,s2) in zip(nstr(mypi,60),nstr(pi,60)):
    print(s1 if s1==s2 else s1+'\u0332',end='')
print()

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

3.14159245358979323846464338327950278419716939938730582097494

и, на самом деле, такого количества слагаемых уже достаточно, чтобы найти не только первые 6, но и первые 60 знаков π


продолжение ниже

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

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

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

просуммируем первые N членов… насколько мы далеки от правильного ответа? сумма оставшихся членов примерно равна интегралу 1/x³ от N до ∞, т.е. ~1/2N²

то есть даже для N=100500 мы получаем всего 10 правильных цифр, а до 100 правильных цифр дойти практически невозможно (требуется ~10^{50} слагаемых)…


…но подождите! «ваш прогноз выполняется с вероятностью 40%? так говорите все наоборот, и будет выполняться с вероятностью 60%» — раз мы знаем, насколько наш ответ далек от правильного, так давайте соответствующим образом скорректируем ответ!

и действительно, для N=100500 сумма первых N членов +1/2N² дает на 5 больше правильных цифр


можно ли оценить сумму (в данном случае, хвост ряда) точнее, чем просто интегралом? оценка «интеграл ≈ сумма» возникает при приближении площади под графиком площадью столбиков — а ясно, что лучше заменить прямоугольники трапециями

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

на самом деле, +1/2N²-1/2N³ — это начало серии поправок, а возникающие коэффициенты — это (примерно) числа Бернулли… писал про этот сюжет (формулу Эйлера-Маклорена и числа Б.) в статье «Алгебра, геометрия и анализ сумм степеней последовательных чисел», https://www.mathnet.ru/rus/mp882

***

посмотрим, как это работает на практике

возьмем, скажем, 57 поправочных членов… а сумму не обязательно даже считать до N=100500, ограничимся N=500 — и уже этого хватает для 100 правильных цифр после запятой!

from mpmath import *
mp.dps = 102

from textwrap import wrap

def bernoulli_upto(m):
    ans = [mpf(1)]
    for k in range(1,m+1):
        b = mpf(0)
        for i in range(k):
            b -= ans[i]*binomial(k+1,i) / (k+1)
        ans.append(chop(b,tol=1e-10))
    return ans

def improve(pre,N,m):
    bernoulli = bernoulli_upto(m-1)
    for k in range(1,len(bernoulli)+1):
        pre += (bernoulli[k-1]*k) / (2*N**(k+1))
    return pre

def partial_sum(N):
    ans = mpf(0)
    for k in range(1,N+1):
        ans += mpf(1) / k**3
    return ans

def mprint(num,prec):
    intp, frcp = nstr(num,prec).split(".")
    print(f"{intp},",*wrap(frcp,width=5))

#N = 100500
#ans = partial_sum(N)
#mprint(ans,21) # всего 10 правильных цифр

N = 500
ans = partial_sum(N)
mprint(improve(ans,N,57),101) # все 100 цифр правильные!

***

так что, может вообще зафиксировать N, а только увеличивать количество поправочных членов (аргумент m функции improve)? степени N экспоненциально убывают, ну какие-то коэффициенты стоят при них…

проблема в том, что если начало последовательности Бернулли выглядит безобидно (1, -1/2, 1/6, -1/30…), то дальше эти числа начинают быстро — почти как факториалы! — расти (и для m=57 появляются коэффициенты ~10^{30})

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

***

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

про это поведение ряда Лейбница — взял из книжки «Mathematics by Experiment» (Borwein, Bailey) — и книжку вообще, и конкретно это место в ней мне показал Сережа Маркелов

а сама математика восходит к Эйлеру — которому очень хотелось посчитать 20 знаков суммы обратных квадратов, не складывая 10^{20} чисел

(после этого Эйлер нашел и точное значение этой суммы… а также ζ(2n), ζ(-k), и даже функциональное уравнение для ζ(x) — видимо поэтому ζ(x) теперь называют дзета-функцией… Римана)

сохраню небольшое видео с введением в Lean:

https://youtu.be/I2zaPoj3G50

мб и что-то такое надо попробовать включить в компьютерные развлечение, не только питон

Рассказывал пару раз на летних школах про количества решений разных уравнений mod p. Вот, например, на такую таблицу можно попробовать посмотреть и поискать какие-то закономерности.

Генерировал ответы в таком духе:

from prettytable import PrettyTable

def isPrime(n):
    if n == 2 or n == 3: return True
    if n%2 == 0 or n<2: return False
    for i in range(3, int(n**0.5)+1, 2):
        if n%i == 0: return False    
    return True

t = PrettyTable(['p','#{x²=-1}','#{x²=2}','#{x²+y²=1}',\
  '#{y²=x³+x²}','#{y²=x³-x}'])
for p in range(3,50+1,2):
 counts = [0]*5
 if isPrime(p):
            for x in range(p):
                if (x*x+1)%p == 0: counts[0] += 1
                if (x*x-2)%p == 0: counts[1] += 1
                for y in range(p):
                    if (x*x+y*y-1)%p == 0: counts[2] += 1
                    if (x*x*x+x*x-y*y)%p == 0: counts[3] += 1
                    if (x*x*x-x-y*y)%p == 0: counts[4] += 1
            t.add_row([p,*counts])
 else:
  t.add_row([p,'','','','',''])
t.align="r"
print(t)

Пример того, на что еще интересно смотреть (и при этом есть достаточно доступный ответ) — количества точек на сферах mod p, начиная с x²+y²+z²=1 mod p.

Если хочется менее расплывчатых вопросов — то вот такие, например, задачи давал: https://dev.mccme.ru/~merzon/leto2022/p-count-ex.pdf

Если, наоборот, хочется возвышенных разговоров — то вот (для школьников, не для экспертов) есть рассказа А.Кузнецова про теорему Ферма, гипотезу Римана и алгебраическую геометрию: https://youtu.be/ElFNGh4WF6Y

в конце поста https://www.group-telegram.com/compmathweekly.com/6 был вопрос без ответа, можно ли так получить асимптотически правильную (пусть с точностью до константы) оценку на факториал

а сегодня коллега Устинов всё объяснил

// дальше техническое математическое, простите

напомню контекст: e^n = …+n^n/n!+…, поэтому e^n>n^n/n! (мы оценили сумму положительных чисел самым большим слагаемым), т.е. n!>(n/e)^n

так вот, если n^n/n!=S(n), то следующий член равен S(n)/(1+1/n), т.е. больше, чем S(n)e^{-1/n}; следующий равен S(n)/(1+1/n)(1+2/n), т.е. больше, чем S(n)e^{-1/n}e^{-2/n} и т.д. — если отойти на m шагов, то получится хотя бы S(n)e^{-m(m+1)/2n}

в частности, если отойти на (как и было предсказано :) на √n шагов, то получится хотя бы с⋅S(n) (для некоторого с, не зависящего от n) — и это дает оценку e^n > √n с S(n), т.е. n! > с√n (n/e)^n

можно доказать и оценку с другой стороны:
(1+1/n)(1+2/n)…(1+m/n)>1+1/n+…+m/n>1+m²/2n, поэтому если отойти на √n шагов, то слагаемое уменьшится хотя бы в полтора раза и т.д. — получается, что n! < C√n (n/e)^n

можно вытащить и полноценную формулу Стирлинга с использованием формулы Валлиса (оценивая центральный биномиальный коэффициент)

пусть мы хотим разбить геометрическую прогрессию 1,q,…,q^{N-1} на 4 равные части E, F, G, H

если E(q)=F(q), то многочлен F-E делится на минимальный многочлен числа q — и то же верно для G-E, H-E

но тогда делится на этот минимальнй многочлен и (F-E)+(G-E)+(H-E)=(E+F+G+H)-4E=1+q+…+q^{N-1}-4E

это условие совершенно не достаточное, только необходимое, но все же попробуем так протестировать гипотетическую долю для котика E=1:

import sympy
q = sympy.Symbol('q')
E = 1
for N in range(3,40):
    print(N,sympy.factor(-4*E+(1-q**N)/(1-q)))

видим, что нужный многочлен получается приводимым в обозримом диапазоне только дважды: для N=4 (этот случай явно не подходит) и для N=14, где он делится на q³+q²-1

теперь можно попробовать взять N=14, доли двух котиков 1 и q²+q³, а доли оставшихся двух котиков поискать прямолинейным перебором (раздаем оставшиеся куски всевозможными способами между G и H и проверяем, делятся ли G-1 и H-1 на q³+q²-1)

про это программу положу в комментарии

***

можно попробовать так же подойти к 5 котикам, но там для E=1 в обозримом диапазоне кандидатов на решение не видно!

можно пытаться перебирать разные E, но тут, наоборот, начинает лезть много паразитных решений (напр., 1+q^7, q+q^4, q^2+q^3, q^5+q^6 все равны при q=-1, но мясо это резать не помогает )

в общем, если сможете найти пример для деления на 5 частей, то пишите

***

впечатления от sympy… неоднозначные

например, заменяю sympy.factor на sympy.factor_list (чтобы сразу получить множители списком) — и сразу все падает с неясными ошибками… и так все время

нравится сюжет Конвея про аналогию между играми и числами

например, игры (скажем, в которых роли противников симметричны, а проигрывает тот, кто не может сделать ход) можно складывать: в G+H играют на двух столах, на одном столе позиция в игре G, на другом — в игре H, каждый раз можно выбрать один из столов и сделать за ним ход

если в H выигрывает второй игрок, то результат у G+H такой же как и в G — это мотивирует объявить все выигрышные для второго игрока игры нулевыми

а вот игры, в которых выигрывает первый, бывают очень разными

если «ним-число» *n — это глуповатая игра «есть кучка из n камней, за ход можно взять любое количество камней из кучки», то *0 действительно нулевая игра, а все остальные *n — различные… и ненулевые )

и игра в четыре кучки камней *1+*3+*5+*7 уже не очень простая (не все персонажи фильма L'Année dernière à Marienbad справились), чтобы научиться в нее играть, хорошо бы изучить таблицу операций с ним-числами

вот такой, например, листок про это: https://dev.mccme.ru/~merzon/v14/pscache/5d-nim.pdf

написал код, который выписывает таблицы сложения и умножения для ним-чисел

def mex(N,arr):
    for a in range(N):
        if (a not in arr):
            return a
    return None

N = 2**(2**2)

t_sum = [list(range(N))]
for m in range(1,N):
    newline = []
    for i in range(N):
        # *m+*i = mex{*j+*i,*m+*i'|j<m,i'<i}
        arr = [line[i] for line in t_sum] + newline
        newline.append(mex(N,arr))
    t_sum.append(newline)
print(*t_sum,sep="\n")

t_mul = [[0]*N]
for m in range(1,N):
    newline = []
    for i in range(N):
        # *m.*i = mex{*j.(*i+*i')+*m.*i'|j<m,i'<i}
        arr = []
        for i1,mi1 in enumerate(newline):
            ii1 = t_sum[i][i1]
            for line in t_mul:
                jii1 = line[ii1] #*j.(*i+*i')
                arr.append(t_sum[jii1][mi1])
        newline.append(mex(N,arr))
    t_mul.append(newline)
print()
print(*t_mul,sep="\n")

можно заметить, а потом и доказать, что ним-сложение — это, на самом деле, просто побитовое сложение

а вот для ним-умножения настолько простого описания, кажется, нет

( определение — можно прочитать в https://en.wikipedia.org/wiki/Nimber#Multiplication )

но операция оч. хорошая — в частности, ним-числа, меньшие *(2^(2^k)), образуют конечное поле

(мини-комментарий к предыдущему)

как расширить поле из 2 элементов до поля из 2^2=4 элементов? достаточно добавить элемент A, удовлетворяющий соотношению
A²=A+1

дальше получить поле из 2^(2^2) элементов можно расширением B²=B+A; потом можно взять расширение C²=C+B и т.д.

именно такую последовательность расширений и реализует ним-умножение (A=*(2^(2^0)), B=*(2^(2^1)), C=*(2^(2^2))…)

(для алгебраического замыкания поля из 2 элементов всего это, конечно, маловато — но, говорят, можно взять все ним-числа до какого-то подходящего ординала, тогда действительно будет алгебраически замкнуто… ну и «все» ординалы как ним-числа образуют алгебраически замкнутое «поле»)