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

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

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

--Г.М.

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/

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