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

Г.Б.Шабат попросил продемонстрировать это на примере уравнения x^3-8x^2+5x+1=0 (дискриминант 7^4) — оказалось, как нередко бывает, что между «в принципе понимаю как делать» и «могу сделать на практике» есть некоторый зазор (но постепеннно все получилось)

1.
легко попросить sympy решить кубическое уравнение:

from sympy import *

x = Symbol('x')
P = x**3 - 8*x**2 + 5*x + 1
theta0, theta1, theta2 = solve(P,x)

но после print(theta0) мы увидим длинную формулу с двумя кубическими корнями из комплексных чисел, прямо из нее непонятно как извлечь ответ с косинусами

2.
чтобы решать кубическое уравнение, удобно сделать преобразование Фурье и смотреть не на корни, а на такие их линейные комбинации («резольвенты Лагранжа»):

omega = exp(2*pi*I/3)
lambda0 = theta0 + theta1 + theta2
lambda1 = theta0 + omega*theta1 + omega**2*theta2
lambda2 = theta0 + omega**2*theta1 + omega*theta2
чем это полезно? тем, что когда группа Галуа переставляет тэты по циклу — lambda1 и lambda2 домножаются на степени ω, соотвественно их кубы неподвижны, т.е. лежат в Q[ω]

собственно, это в формуле для корней кубического уравнения мы и видим: theta0=(lambda0+lambda1+lambda2)/3, и первое слагаемое — коэффициент нашего уравнения (в данном случае, 8), а вторые два — кубические корни каких-то крокодилов, лежащих в Q[ω]

3.
как найти нужного крокодила? print(lambda1**3) дает невразумительное выражение на несколько строк, а хотелось бы видеть конкретный ответ a+bω с рациональными (на самом деле, даже целыми) a и b… спросим sympy получше:

print(minimal_polynomial(lambda1,x,domain=QQ))

x**6 - 637*x**3 + 117649
число 117649 может пугать (если не сразу заметили, что это 7^6) — но в любом случае, мы получили просто квадратное уравнение (на куб lambda1), его можно решить и получить -7²(1+3ω)²

(можно на всякий случай проверить: print(N(lambda1**3+7**2*(1+3*omega)**2)) дает более-менее ноль)

если читали https://t.me/compmathweekly/23 — можно подумать, как теперь получить конкретный ответ )

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

бонусная картинка в продолжение предыдущего: если простое число представимо в виде суммы а) 2; б) 3; в) 4; г) 8 квадратов, то сколько таких представлений

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

(там еще добавлены для мастштаба две линии типа с√n — но на самом деле максимумы растут примерно как √n log log n, насколько понимаю)

1.
понравился прием с reshape — покажу на примере: пускай мы интересуемся, скажем, представлением чисел в виде суммы квадратов


import numpy as np

def sum_squares(N,k):
counts1 = np.zeros(N*N+1,dtype=int)
counts1[0] = 1
np.add.at(counts1,[n*n for n in range(1,N+1)],2)
counts = counts1.copy()
for _ in range(k-1):
counts = np.convolve(counts,counts1)
return counts[:N*N+1]
(такой подсчет сумм сверткой — это примерно как возводить производящую функцию квадратов в степени, только без “развлечений” с символическими вычислениями)

тогда можно печатать ответы по сколько-то в строку в духе

print(np.reshape(sum_squares(20,3)[:240],(-1,8)))
(число 8 можно подобрать, желая чтобы расположение нулей выглядело не так хаотично — и возникающую гипотезу дальше легко доказать)


2.
для сумм 4, скажем, квадратов ничего из reshape вроде бы не видно… зато помогает посмотреть на значения для простых чисел — на бумажке это немного утомительно, а компьютерно легко (все числа там получаются кратными 8, так что сразу на 8 разделим):


counts = sum_squares(N:=10,4)
for p in np.nonzero(is_prime(N*N))[0]:
print(f"{p}: {counts[p]//8}")

какая возникает гипотеза? )

а что будет для p², p³? (спойлер: ответ равен количеству одномерных подпростравнств в 2-, 3-, 4-мерных пространствах над Z/p — это приятно понимать при помощи кватернионов)

(upd)
заодно посчитал количество представлений простых в виде сумму 8 квадратов — и тут легко угадал и неизвестный мне заранее ответ: 16(p³+1) — мб надо подумать, как это доказать

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

пусть A — верхнетреугольная матрица, с которой мы начинали (единицы стоят там, где j делит i), тогда M(n) — это просто значение первой переменной в решении системы Ax = (1…1), посчитанное при помощи правила Крамера

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

A^t — это матрица суммирования по делителям; обратная операция к суммированию по делителям хорошо известна — это свертка с функцией Мёбиуса (μ=(-1)^s для произведения s различных простых, 0 иначе)

отсюда получаем, что M(n) — это сумма μ(k) по k от 1 до n

в частности, ясно, что |M(n)| не превосходит n (т.к. |μ|⩽1); кстати, оценка |M(n)|=o(n) эквивалентна теореме о распределении простых чисел (количество простых до n ~ n/log(n)) — в общем, связь этих определителей с аналитической теорией чисел перестает казаться такой уж загадочной

но тут самое время сделать паузу в разговорах и это вычисление M(n) реализовать:


import numpy as np
import matplotlib.pyplot as plt

def primes(N):
is_prime = np.ones(N+1)
is_prime[:2] = 0
for n in range(int(np.sqrt(N))+1):
if is_prime[n]:
is_prime[n*n::n] = 0
return np.nonzero(is_prime)[0]

N = 3*10**6

moebius = np.ones(N+1)
for p in primes(N):
moebius[::p] *= -1
moebius[::p**2] = 0
mertens = np.cumsum(moebius)
mertmax = np.maximum.accumulate(np.abs(mertens[1:]))

ax = plt.axes()
ax.plot(range(1,N+1),mertens[1:])
ax.plot(range(1,N+1),np.sqrt(range(1,N+1)),-np.sqrt(range(1,N+1)))
ax.plot(range(1,N+1),mertmax,-mertmax)
plt.xlim(0,None)
plt.show()

функция Мёбиуса принимает значения ±1 довольно хаотичным образом, поэтому можно ожидать чего-то в духе закона повторного логарифма (который выполнялся бы, если бы mu действительно генерировалась подкидываниями монетки), что максимум |M(n)| мог бы расти примерно как √n × (log log n)^(1/2)… есть гипотеза, что в реальности он растет примерно как √n × (log log log n)^(5/4) (а гипотеза Римана эквивалентна оценке √n O(n^ε) для всех ε>0)

(ситуация напоминает обсуждавшуюся в https://t.me/compmathweekly/18 — с математической точки зрения кратные логарифмы стремятся к бесконечности, но увидеть это в компьютерном эксперименте практически невозможно, уже двойной логарифм везде выглядит как ограниченная функция)

и доказано, что для каких-то n функция Мёртенса M(n) принимает значение, большее по модулю чем √n… но оценка сверху на первое такое n — что-то типа 10^(10^60) (!), а для всех реально вычисленных значений |M(n)|/√n < 0.6

рассмотрим матрицу n×n, в первой строке которой одни единицы, во второй — 0101…, в третьей — 001001… и т.д.

можете посчитать ее определитель? ну да, она верхнетреугольная, определитель 1

теперь еще в каждой строке первый элемент заменим на 1 — чему равен определитель M(n) получившейся матрицы? ну или по крайней мере можете ли доказать на него оценку, что по модулю он меньше √n?

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

были поставлены разные эксперименты, и для всех вычисленных значений |M(n)| не превосходит 0,6√n (вот на картинке фрагмент экспериментальных данных)

via И.Тайманов; продолжение следует

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

работа над ошибками

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


import flint

def count_tilings(n,m,ext):
ans = flint.fmpz_mat(2**n,1)
ans[0,0] = 1
for _ in range(m):
ans = ext * ans
return ans[0,0]

def newext(n,ext):
ans = flint.fmpz_mat(2**n,2**n)
for mask in range(2**n):
for perp in range(2**n):
if mask&perp == 0:
if (mask+perp)%2 == 1:
ans[mask,perp] = ext[n-1][mask>>1,perp>>1]
elif (mask+perp)%4 == 0:
ans[mask,perp] = ext[n-2][mask>>2,perp>>2]
return ans

N = 5
ext = [flint.fmpz_mat(1,1,[1]),flint.fmpz_mat(2,2,[0,1,1,0])]
for n in range(2,2*N+1):
ext.append(newext(n,ext))
if n%2 == 0:
ans = count_tilings(n,n,ext[n])
print(n,ans)

12988816 разбиений квадрата 8×8 теперь подсчитываются за пару сотых секунды



теперь легко проверить, что разные делимости есть и для разбиений прямоугольников на доминошки — например, для 6×2 разбиений 13 => для (7k+6)×(3l+2) количество разбиений делится на 13 (и действительно, для 13×20 получается 5422065916132126528319352874496 разбиений, это число делится на 13)

количество разбиений прямоугольника 2×N — это просто число Фибоначчи, и что F_N | F_{2N+1} легко доказать комбинаторно (посмотреть как там покрыт доминошками средний столбец)

а как доказать комбинаторно аналогичное утверждения для разбиений прямоугольника, say, 6×N — я не знаю (само утверждение тоже прочитал у Проппа)

***

еще переписал генерацию случайного разбиения на доминошки ацтекского брильянта — через domino shuffling (aka urban renewal)

полный код не влезает в лимит телеграма, так что положу только основную функцию, которая делает из случайного брильянта порядка size случайный брильянт порядка size+1… point, наверное, в том, что получилось не сложнее, чем наивный код с флипами (но, конечно, флипы работают для разных областей, а не только для брильянтов)


def extend(field,size):
newfield = [ [' ']*(2*size+2) for _ in range(2*size+2) ]
for i in range(2*size+2):
for j in range(2*size+2):
if inside(i,j,size+1) and newfield[i][j]==" ":
if inside(i,j-1,size) and field[i][j-1] in "Uu" and (not inside(i-1,j-1,size) or not field[i-1][j-1] in "Dd"):
newfield[i][j] = field[i][j-1]
elif inside(i-2,j-1,size) and field[i-2][j-1] in "Dd" and (not inside(i-1,j-1,size) or not field[i-1][j-1] in "Uu"):
newfield[i][j] = field[i-2][j-1]
elif inside(i-1,j,size) and field[i-1][j] in "Ll" and (not inside(i-1,j-1,size) or not field[i-1][j-1] in "Rr"):
newfield[i][j] = field[i-1][j]
elif inside(i-1,j-2,size) and field[i-1][j-2] in "Rr" and (not inside(i-1,j-1,size) or not field[i-1][j-1] in "Ll"):
newfield[i][j] = field[i-1][j-2]
else:
bit = getrandbits(1)
newfield[i][j], newfield[i][j+1] = ("U", "u") if bit else ("L", "R")
newfield[i+1][j], newfield[i+1][j+1] = ("D", "d") if bit else ("l", "r")
return newfield

коллега Бахарев спрашивал чуть раньше, получается ли что-то интересное, если количества замощений доски N×N (или что-то подобное) рассматривать по разным модулям

обычно ничего особо не видно, но прочитал недавно у Проппа, что поучительно рассматривать разбиения квадрата 2N×2N на доминошки и квадратики 2×2 по модулю степеней двойки

и действительно, компутер быстро посчитал мне начало последовательности
1
3
165
205879
5811552169
3676802302990923
51863128477992833745597 16264377685372072288158612804415

выглядит устрашающе, но по модулю 8, скажем, это просто
1, 3, 5, 7, 1, 3, 5, 7…

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

(технически модифицировал чуть-чуть как раз код для подсчета разбиений на доминошки динамикой по профилю из https://t.me/compmathweekly/17 — и немножко лениво сделал, мб из-за этого не очень быстро работает)

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

один из источников — статья Проппа «Some 2-adic conjectures concerning polyomino tilings of Aztec diamonds» (https://arxiv.org/abs/2204.00158)

надеюсь, что «гостевые посты» коллег здесь — тоже станут традицией

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

Как известно, сфера с g ручками получается из 4g-угольника попарной склейкой сторон. Возьмем (при g >= 2) 4g-угольник в плоскости Лобачевского, у которого отождествленные пары сторон равны по длине, а сумма углов равна 2π. Тогда на склеенной из него поверхности остается гиперболическая метрика; обратно, любую гиперболическую поверхность можно разрезать (по геодезическим отрезкам) на 4g-угольник, вершины которого сойдутся в выбранной на поверхности точке.

Если теперь мы рассмотрим универсальное накрытие поверхности плоскостью Лобачевского, то соответствующий 4g-угольник замостит плоскость своими копиями. Так мы можем изображать гиперболические метрики на поверхности (точки в пространстве Тейхмюллера) при помощи замощений.

Я написал небольшую программу, которая рисует фрагмент случайного такого замощения для g=2, т.е. копиями восьмиугольника. Склеиваемые стороны обозначены одинаковым цветом.

Я существенно опирался на архивные материалы аналогичного (но имевшего больший интерактивный функционал) проекта "The Teichmüller Navigator", созданного в 1993 Deva Van Der Werf под руководством David Ben-Zvi и Paul Burchard - каким-то чудом страница проекта сохранилась до сих пор: http://www.geom.uiuc.edu/apps/teich-nav/report/report.html.

ссылка на github с исходным кодом и инструкциями: https://github.com/mntroshkin/teichmuller-explorer.

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

как расширить поле из 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 элементов всего это, конечно, маловато — но, говорят, можно взять все ним-числа до какого-то подходящего ординала, тогда действительно будет алгебраически замкнуто… ну и «все» ординалы как ним-числа образуют алгебраически замкнутое «поле»)

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

например, игры (скажем, в которых роли противников симметричны, а проигрывает тот, кто не может сделать ход) можно складывать: в 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 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 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)), образуют конечное поле

пусть мы хотим разбить геометрическую прогрессию 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 (чтобы сразу получить множители списком) — и сразу все падает с неясными ошибками… и так все время

на Московской математической олимпиаде была сегодня такая задача
«У хозяйки есть кусок мяса, которым она хочет накормить трёх котиков. Раз в несколько секунд хозяйка отрезает кусочек мяса и скармливает его одному из котиков на свой выбор, причём каждый кусочек должен составлять одну и ту же долю куска, от которого его отрезают. Через некоторое время хозяйка убирает остаток мяса в холодильник. Может ли она скормить котикам поровну мяса?» (А.Кушнир)

другими словами, можно ли разбить члены какой-нибудь конечной геометрической прогрессии с 0
любая хозяйка, которая умеет резать мясо в золотом сечении, решит задачу для 2 котиков; для 3 котиков попробуйте решить на бумажке; для 4 котиков покажу попозже (пусть пока будет возможность подумать), как воспользоваться компьютером; какой ответ для 5 котиков — хотел бы знать

в конце поста https://t.me/compmathweekly/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

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

Рассказывал пару раз на летних школах про количества решений разных уравнений 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

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

https://youtu.be/I2zaPoj3G50

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

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

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

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

просуммируем первые 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) теперь называют дзета-функцией… Римана)

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


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 знаков π


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

в дополнение к предыдущему — количества 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)

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

сколькими способами число 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) по разным модулям

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

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

это может показаться смутно гармоничным (если в плоском разбиении на клетке (i,j) лежит кубик, то и на клетках (i',j) и (i,j') при всех i'
это тоже магия соответствия RSK


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

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


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

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

***

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