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

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

Как известно, сфера с 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.

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

коллега Бахарев спрашивал чуть раньше, получается ли что-то интересное, если количества замощений доски 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://www.group-telegram.com/compmathweekly.com/17 — и немножко лениво сделал, мб из-за этого не очень быстро работает)

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

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

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

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

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

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

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

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

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

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

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

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

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

пусть 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://www.group-telegram.com/compmathweekly.com/18 — с математической точки зрения кратные логарифмы стремятся к бесконечности, но увидеть это в компьютерном эксперименте практически невозможно, уже двойной логарифм везде выглядит как ограниченная функция)

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