хочу научиться компьютерным экспериментам в математике

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

начиная с чего-то совсем простого, а там посмотрим, как пойдет

--Г.М.

def prime_factors(num):
    factor = 2
    while num > 1:
        if num % factor == 0:
            yield factor
            num = num // factor
        else:
            factor += 1

for k in [2,3,5,7,11,13,17,19,23,29]:
    print("2^",k,"-1",end="=",sep='')
    print(*prime_factors(2**k-1),sep='x')

в статью в Квантике про числа Мерсенна (вошла в №1 за 2025 год) захотелось включить короткую программу, которая демонстрирует, что числа вида 2^p-1 не всегда являются простыми

2^2-1=3
2^3-1=7
2^5-1=31
2^7-1=127
2^11-1=23x89
2^13-1=8191
2^17-1=131071
2^19-1=524287
2^23-1=47x178481
2^29-1=233x1103x2089

e^x = 1+x+x²/2+x³/3!+…

пусть x какое-нибудь конкретное число (на картинке x=20) — какие слагаемые дают основной вклад в сумму?

нетрудно понять, что для e^n самое большое слагаемое — это n^n/n!

дальше слагаемые уменьшаются как минимум в геометрической прогрессии и можно оценить
n^n/n! < e^n < (2n+1) n^n/n!
что дает простое доказательство слабой формы формулы Стирлинга:
(n/e)^n < n! < Cn (n/e)^n

прикольно, что это не так далеко от правильной асимптотики n! ~ C√n (n/e)^n

наверное оценку можно усилить, если понять какой ширины «горб» на картинке (в духе цпт и т.п.), но можно ли сделать это не опираясь на Стирлинга? не продумывал еще

===

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

файл в комментариях

Не планировал еще одного поста так скоро, но пришел Костя Кноп и рассказал забавное:

Как быстро посчитать много знаков числа π, если есть калькулятор с тригонометрическими функциями (но без обратных тригонометрических)? Предлагается такой рецепт: начинаем с грубого приближения π≈3, а дальше делаем пару-тройку итераций x→x+sin(x).

Можно было бы не писать никакой программы — но раз уж тут канал про компьютерные эксперименты:

from mpmath import *
mp.dps = 200

def iterate(approx):
    return approx+sin(approx)

mypi = 3
for k in range(4):
    mypi = iterate(mypi)
    diff = pi-mypi
    print(k+1,":",nstr(mypi,50),
          "diff:",nstr(diff,2))

Кто запустил код — может видеть, что через 4 итерации ошибка уже порядка 10^{-100} (!)

Ну и объяснить это не сложно: за итерацию мы заменяем π+t на π+t-sin(t), а t-sin(t)≈t³/6 — вот и увеличивается за каждую итерацию количество правильных цифр примерно втрое.

Еще можно заметить, что мы недалеко ушли от метода Ньютона (буквально метод Н. для sin(x) — это x→x-tg(x), но рядом с π это примерно то же).

Тему вычисления знаков π и т.п. наверное еще продолжим.

===

Я ничего не понимаю в библиотеках питона, а mpmath — то что сходу нашлось в гугле по моему запросу (чтобы можно было много цифр синуса считать и т.п.).

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

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

===

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

про прямоугольники 2×N и числа Фибоначчи и так все знают

а вот на картинке ацтекский бриллиант порядка 4 — сколькими способами его можно разбить на доминошки? а бриллиант порядка N?

тут общий ответ можно угадать по ответам для небольших N, но перебирать руками тяжеловато (для N=4 вариантов уже больше тысячи), компьютерный эксперимент действительно полезен

про разные способы этот ответ доказывать есть интересная брошюра Е.Смирнова «Три взгляда на ацтекский бриллиант», https://mccme.ru/free-books/dubna/smirnov-aztec.pdf

***

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

все же перебрал на компутере и замощения квадрата 8×8 — в качестве ответа получилось название лекции С.Смирнова https://www.mathnet.ru/present18250

на картинке — генерация случайного разбиения ацтекского брильянта порядка 40 на доминошки (10 тыс. попыток флипов, 100 тыс., 1 млн., много млн¹)

начал с «большого числа» 10000 итераций, а когда дошел до 1000000 — стало казаться, что ничего не получится… а на самом деле это уже почти успех, важно только не опускать руки

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

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

ну и это еще 40-брильянт не такой уж, на самом деле, большой — в комментарии положу картинку для порядка 100… когда она догенерируется

картинок брильянтов можно найти сколько угодно в интернетах, в «Мат. байках» Витя Клепцын объясняет качественно почему так получается² — но всё равно когда сам запускаю программу и реально этот полярный круг возникает, ощущаю чудо

¹ в середине века у Роллс-Ройсов в техпаспорте в графе «мощность двигателя» было написано просто «достаточная» («adequate»)

² см. https://news.1rj.ru/str/mathtabletalks/733 и далее

===

в следующий раз, наверное, будет что-то low tech — мб с экселем вместо питона

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

на картинке — компьютерный эксперимент на эту тему: в экселе легко¹ сгенерировать таблицу для сумм двух квадратов

если N=X²+Y², то 2N=(X+Y)²+(X-Y)², поэтому интересна прежде всего представимость нечетных чисел²

вот под основной таблицей посчитано, какие из нечетных чисел до 20 представимы в виде суммы двух квадратов (встречаются в основной таблице), а какие нет

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

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

¹ если нужны подробности, как это сделать — см. комментарии

² вместо алгебры можно думать про геометрию (а вместо экселя можно, конечно, пользоваться питоном) — см., например, статью про косые квадраты в https://kvantik.com/issue/pdf/2021-07.pdf

сохраню ссылку на семинар по экспериментальной математике, который 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

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