Компьютерная математика Weekly

краткий апдейт на тему t.me/compmathweekly/141

just for fun на каникулах: purplesyringa.moe/blog/log-is-non-monotonic-in-php-and-lua/ — реализация логарифма иногда оказывается немонотонной функцией (и не по той причине, о которой вы, вероятно, подумали)

история про Rowland'а и Sinkhorn limit немного повисла в воздухе — вернемся ненадолго

матрица 2×2, у которой суммы по строкам и столбцам единицы, имеет вид [x 1-x\\ 1-x x]

чтобы понять, к какой из матриц такого вида мы сойдемся, нужно еще найти инвариант процесса — ну вот он такой соверешнно в духе мат. кружка: произведение чисел на черных клетках делить на произведение чисел на белых, ad/bc — получается уравнение на x, из которого x = √ad/(√ad+√bc) — вот эти квадратные корни и были видны в эксперименте

конечно аналогичные инварианты можно писать для 2×2 подматриц матрицы 3×3 — и в итоге на элементы предельной матрицы получаются алгебраические уравнения степени 6

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

будем переходить от многоугольника к новому многоугольнику с вершинами в серединах сторон исходного («отображение Бюффона»)

если итерировать это отображение, то все вершины стремятся к одной точке… но что при этом происходит с формой N-угольника?

оказывается, он становится аффинно правильным (это, кстати, достаточно просто доказать); в приземленных терминах: после нескольких итераций он практически становится «параллельником» — (практически) параллельны те же стороны и диагонали, что у правильного N-угольника

(для 4-угольников уже после первой итерации получается параллелограмм… ну вот Бюффона можно считать таким обобщением Вариньона)

можно посмотреть пример на гифке или попробовать разные многоугольника по ссылке dev.mccme.ru/~merzon/compmath/midpoints.html

***

услышал про это в начале лекции А.П.Веселова на ЛШСМ-2026 — дальше обсуждался трехмерный случай, где все сложнее

upd: появилось, кстати, видео — mathnet.ru/rus/present50983

во время ЛШСМ на компьютерные развлечения не хватает энергии, так что вот пока вместо моего поста — пост Тао:

terrytao.wordpress.com/2026/07/14/visualizing-the-gilbreath-expectation-sequence/

упомянутый в прошлом посте Rowland (относительно) недавно рассказывал, оказывается, на семинаре по экспериментальной математике вот про что

возмьем квадратную матрицу неотрицательных чисел. и будем нормировать строки-столбцы: разделим каждую строку на сумму чисел в ней, потом каждый столбец, потом снова каждую строку… к чему это сойдется («Sinkhorn limit»)?

вот, например, для матрицы [4 1\\ 2 1] можете сообразить, что это за числа получаются? можно посмотреть на приближенные значения:

import numpy as np
A = np.array([[4,1],[2,1]])
for _ in range(5):
    A = A / A.sum(axis=1, keepdims=True)
    A = A / A.sum(axis=0, keepdims=True)
    print(A)

(для матриц произвольного размера жизнь быстро усложняется — и доклад был как раз про разные экспериментальные гипотезы по этому поводу… мб напишу позже какие-то подробности)

рассмотрим последовательность
a(1) = 7
a(n) = a(n-1)+НОД(n, a(n-1))
Rowland доказал, что каждый раз число увеличивается либо на 1, либо на простое число (но появляются ли так все нечетные простые, неизвестно)

( и такая задача предлагалась, как научили в комментариях, на Турнире городов — problems.ru/view_problem_details_new.php?id=64532 )

коллега Медведь поделился забавным родственником этой последовательности, который для разных начальных условий (гипотетически) генерирует всё большие простые числа-близнецы:

https://math.stackexchange.com/q/5142627/152

(по ссылке как раз эксперименты и обсуждения)

рассказывал сегодня про треугольник Серпинского — и хотелось показать, как он возникает в «игре в хаос»

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

а если запустить много кузнечиков, то треугольник Серпинского постепенно появляется на экране

позапускать кузнечиков можно по ссылке dev.mccme.ru/~merzon/compmath/sierpinski.html

возьмем теперь не конкретные многочлены, а что-то более случайное…

ну например, будем брать многочлены степени 200 со старшим коэффициентом 1, а остальные к-ты будем выбирать… ну скажем равномерно из отрезка [-10;10]

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

см. также: mathoverflow.net/q/182412/1556

в качестве интермедии — нарисовал просто картинки разных многочленов на комплексной плоскости

цвет соответствует аргументу значения в данной точке

хорошо видны корни — точки, рядом с которыми сходятся все цвета

это само по себе намекает на одно из доказательств основной теоремы алгебры

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

ну так… с трудом (много жестких конфигураций, всё постоянно застревает где-то не там). у меня после серии доработок какие-то упаковки получается увидеть, но (даже для небольшого числа квадратов) не наилучшие

не хватает каких-то (математических) идей

оставлю просто пару картинок

в квадратную коробку какого наименьшего размера можно положить N единичных кругов?

если N большое, то известно что примерно мы увидим… или если N какое-нибудь круглое, типа 4 или 9… а если N какое-нибудь дурное, типа 11?

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

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

в simulated annealing эти две идеи смешаны: сначала температура высокая и мы делаем достаточно случайные переходы, потом температура понижается и ухудшающие ситуацию переходы становятся всё менее вероятными… и всё магическим образом работает

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

вот, собственно, практически весь код:

def simulated_annealing(N, max_iter, T=0.2, cooling=0.9975):
    centers = np.random.uniform(0, 1, (N, 2))
    best_centers = centers.copy()
    best_R = R = max_radius(centers)
    history = [R]

    for step in range(max_iter):
        i = step%N
        old_pos = centers[i].copy()
        old_dist = point_dist(i, centers)

        step_size = min(T**0.5,0.04)
        centers[i] += np.random.normal(0, step_size, 2)
        centers[i] = np.clip(centers[i], 0.0, 1.0)

        delta = point_dist(i, centers) - old_dist
        if delta > 0 or np.random.random() < np.exp(delta / T):
            R = max_radius(centers)
            if R > best_R:
                best_R = R
                best_centers = centers.copy()
        else:
            centers[i] = old_pos

        T *= cooling
        history.append(R)

    return best_centers, best_R, history

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

на наилучшие известные упаковки можно посмотреть на странице erich-friedman.github.io/packing/ — мб кто-то из читателей сможет что-то из рекордов улучшить ;

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

для точки P=(3,5) на кривой y²=x³-2 считали уже координаты точек nP. можно заметить, что знаменатели иксов все время оказывались точными квадратами… ну вот найдем, чьи именно это квадраты:

from fractions import Fraction
import math

N = 20

x, y = x0, y0 = 3, 5
u = [0]*(N+1)
u[1] = x0.denominator
for n in range(2,N+1):
    k = Fraction(y-y0, x-x0) if (x,y) != (x0,y0) \
            else Fraction(3*x0*x0, 2*y0)
    x = k*k-x0-x 
    y = -(k*(x-x0)+y0)
    b = math.isqrt(b2 := x.denominator)
    assert b**2 == b2
    u[n] = b

получается последовательность 1, 10, 171, 7660, 12660211, 22652313570…

так вот, можно рядом написать рекурренту типа «Сомос-4» (v[n+2]·v[n-2] = c₁·v[n+1]·v[n-1]+c₂·v[n]² — и рекуррентам именно такого вида удовлетворяют также числа Фибоначчи или количества замощений ацтекского брильянта):

v = [0]*(N+1)
v[:5] = [0, 1, 10, 171, -7660]
c1, c2 = v[2]**2, -v[3],
for n in range(5,N+1):
    b = c1*v[n-1]*v[n-3]+c2*v[n-2]**2
    assert b%v[n-4] == 0
    v[n] = b//v[n-4]

    print(f"n = {n:2d}: {u[n] == v[n] or u[n] == -v[n]} ({len(str(u[n])):3d} digits)")

и убедиться, что всё сходится


как введение в последовательности Сомоса в hands-on стиле — был проект на ЛКТГ-2023 и статья А.Устинова в Кванте №№8-9 за 2023 год

несколько пренебрегая принципом «show, don't tell», хотел кратко написать про связи (местами пунктирные) между некоторыми из сюжетов здесь

начнем с конца. для рациональной точки P на эллиптической кривой знаменатель nP растет примерно как c^{n²}

раньше обсуждались замощения доминошками области на плоскости… и там часто количество замощений растет с той же асимптотикой, c^{площадь}

например, для обсуждавшегося ацтекского брильянта ответ — 2^{n(n+1)/2}. этот ответ можно «сконденсировать», доказав рекурренту M(n+1)M(n-1)=2M(n)²

бывают разные квадратичные рекурренты в таком духе, в т.ч. упоминавшиеся здесь мельком знаменитые последовательности Сомоса… и, скажем, Сомос-4, действительно, кодирует сложение на эллиптической кривой

у этого всего есть игрушечные версии: можно мостить не по-настоящему двумерную фигуру, а более-менее одномерную — прямоугольник 2×N (или 3×N и т.п. — такого рода вещи где-то в начале обсуждались), тогда ответы получаются типа Фибоначчи, которые удовлетворяют [не только квадратичным, но и] линейным рекуррентам, имеют более простую асимптотику c^n

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

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

разные более конкретные вещи тоже можно пытаться переносить: скажем, F_n | F_{nm} — и вот для последовательности знаменателей nP (скажем, сгенерированных кодом из предыдущего поста конкретно) верно буквально то же… и т.п.

незаконченное обсуждение арифметико-геометрического среднего тоже связано со сложением на кубике, AGM реализует «эллиптический логарифм» (т.е., наоборот, позволяет по точке xP найти x… вещественное или даже комплексное)

но пока step into the elliptic realm не выходит, только трогаю пальцами холодную воду

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

вот на последнее как раз легко посмотреть экспериментально

конечно сложение точек уже реализовано в sage и т.п., но по сути это просто проведение секущей через пару точек на кубике и нахождение третьего пересечения с кубикой — и это легко сделать и руками:

from fractions import Fraction

# y^2 = x^3 - 2
def add(P, Q):
    x1, y1 = P
    x2, y2 = Q

    k = Fraction(y2 - y1, x2 - x1) if P!=Q \
            else Fraction(3*x1*x1, 2*y1)

    x3 = k*k - x1 - x2 # Vieta
    y3 = -(k*(x3 - x1) + y1)

    return (x3, y3)

P = (3, 5)

Q = P
for n in range(2,14):
    Q = add(Q, P)
    print(f"[{n:2d}] {Q[0]}")

экспоненциальный рост прямо визуально виден

выписанные числа на экране образуют параболу, то есть рост ~exp(cn²)


ср. это с тем, что происходит для кривой с похожим уравнением, y²=x³+x², но особой (имеющей самопересечение)… — арифметика, как уже говорилось, помнит про геометрию

пусть кривая (скажем) задана полиномиальным уравнением на x и y с целыми коэффициентами. есть общий принцип в духе «арифметика отражает геометрию», что количество решений что-то помнит про геометрию этой кривой

слова про «количество решений» при этом можно понимать по-разному. можно считать решения по разным простым модулям — t.me/compmathweekly/45 & t.me/compmathweekly/69

но можно и более прямолинейно: считаем рациональные решения (X/Z,Y/Z) для |X|,|Y|,|Z|<N и смотрим на асимптотику по N (чтобы не считать одно и то же много раз, желательно еще запретить X, Y, Z иметь общий делитель)

для степени кривой возможны такие варианты: 1-2, 3, много

1-2) для разминки можно подумать, что будет для линейного уравнения? а для x²+y²=1? а для другой коники?

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

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

3) для эллиптической кривой (гладкой кубики) рост логарифмический, причем в асимптотике виден и ранг — всё растет примерно как c(log N)^(r/2)

много) дальше уже никакой асимптотики нет: по теореме Фальтингса при большей степени [если кривая гладкая] — число рациональных точек вообще конечно

***

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

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

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

мб будет еще продолжение

Черепашка Фурье и черепашка Гаусса

Была такая популярная форма программирования для начинающих «черепашка». У нас будет базовая версия, которая умеет 1) двигаться прямо, оставляя след; 2) поворачиваться (на месте) на указанный угол.

Можно написать что-нибудь типа

import math
import matplotlib.pyplot as plt

x, y, phi = 0, 0, 0

def move(s=1,color='blue'):
    global x, y
    x0, y0 = x, y
    x = x0 + s*math.cos(phi)
    y = y0 + s*math.sin(phi)
    plt.plot([x0,x],[y0,y],color=color)

def turn(s):
    global phi
    phi += s*2*math.pi

for k in range(100):
    move()
    turn(1/100)

plt.gca().set_aspect('equal')
plt.show()

и увидеть на экране окружность. Ну… тоже приятно, но не особо интересно.

Забавно, что содержательные вещи буквально в одном шаге от такого. Например, можно вместо move() написать move(<функция от k>) и смотреть, как черепашка занимается преобразованием Фурье.

Но мы лучше оставим сдвиги одинаковыми, а вот поворот пусть растет линейно со временем. Вот на картинке результат

p = 101
for k in range(p):
    move()
    turn((2*k+1)/p)

Можно задавать вопросы… например, далеко ли черепашка уйдет от начала координат при разных p? Тут есть несколько уровней 1) экспериментальный; 2) эмпирическая прикидка в вероятностном духе для большого p; 3) точный ответ для произвольного p.

Можно, кстати, от 2) пойти чуть-чуть в другую сторону и смотреть для больших p на всю форму кривой (а не только на то, далеко ли мы уйдем). Видим мы, на самом деле, спирали Корню (=кривизна пропорциональна пройденному расстоянию).

***

Напомнил о таком сюжете недавно коллега Гусарев, а еще чуть раньше про это писал Гаусс коллега Клепцын.

в новом Квантике (№5) — статья Олега Мухаметова про сумму цифр степеней числа, почему она становится большой

по этому поводу не мог не поставить мини-эксперимент

будем считать сумму двоичных (для вычислительной эффективности) цифр степеней… ну хотя бы 3 — что можно ожидать увидеть?

ну если так грубо, то n log₂3 цифр числа 3ⁿ довольно случайные, так что единиц среди них должна быть примерно половина

вот какие-то колебания рядом с cn для c = (log₂3)/2 мы и видим на графике

import matplotlib.pyplot as plt
from math import log

ns = range(3_000)
xs = [pow(3,n) for n in ns]
ans = [x.bit_count() for x in xs]

c =  log(3)/(2*log(2))
appr = [n*c for n in ns]

plt.plot(ns,ans)
plt.plot(ns,appr)
plt.title(r'2-digit sums of $3^n$')
plt.tight_layout()
plt.show()