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

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

перебирал еще на компутере кое-что из задач, рассматривавшихся для Матпраздника, но подобностей, естественно, не будет — вместо этого будет подборка разного


1.
в обсуждении представимости чисел в виде суммы двух квадратов Ф.Бахарев поднял тему эффективного нахождения таких представлений — а оказывается, про это была статья в Мат. Просвещении: https://www.mathedu.ru/text/mp_2006_v10/p190/

для алгоритма из статьи нужно, правда, сначала найти корень из -1 mod p — это Федя так учит делать:

проще всего найти квадратичный невычет и возвести его быстрым возведением в степень (p-1)/2. невычет можно либо случайно искать (вероятность 1/2), либо просто перебирать первые простые и считать символ Лежандра. если гипотеза Римана верна, то найдется очень быстро. но на практике вообще среди первых нескольких простых обычно есть

а Р.Гусарев поделился статьей про быстрое нахождения представления любого числа в виде суммы 4 квадратов: https://campus.lakeforest.edu/trevino/finding4squares.pdf (и там не обходится, естественно, без кватернионов)


2.
в обсуждении суммы обратных простых М.Трошкин рассказал, что если складывать обратные только к простым числам-близнецам, то сумма сходится (хотя и очень медленно)

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

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


3.
Л.Петров у себя в канале делится демонстрацией того, как себя ведут собственные значения случайных матриц: https://lpetrov.cc/simulations/2025-01-28-goe/

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


наконец, у канала (и прикрепленного к нему чата) появился юзерпик — спасибо Тане Кор. за картинки

Channel photo updated

0.
здесь уже обсуждалось, что если функция достаточно хорошая, то уравнение f(x)=0 можно быстро приближенно решать при помощи метода Ньютона

если x — приближенное значение корня, то рядом ним график функции недалеко ушел от касательной, поэтому в качестве следующего приближения можно взять пересечения касательной с нулем, т.е. x→x-f(x)/f'(x)

и в некоторой окрестности корня метод Ньютона, если производная в этом корне не равна 0, сходится очень быстро: за итерацию погрешность ~возводится в квадрат

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

на картинках выше показано, как работает метод Ньютона для нахождения корня уравнения z⁴-1=0 [на комплексной плоскости]

(какой-нибудь код будет в комментариях)

1.
на левой картинке — сколько итераций надо сделать, чтобы подойти к корню ближе чем на 10^{-6}; видно, что обычно действительно достаточно пары итераций — но иногда даже для вроде бы небольших чисел требуются многие десятки итераций

возникает вопрос, так ли хорош метод Ньютона

Грубо говоря, Немировским было показано, что метод Ньютона локально является оптимальным. Под «локально» имеется в виду, что точку старта метода нужно брать достаточно близко к решению. А что будет, если точку старта (начальное приближение) не удастся подобрать достаточно близко к решению? Ведь мы не знаем, где решение. Собственно, в этом же и задача — найти решение! Какой метод будет оптимальным без таких оговорок?

цитата — из статьи А.Гасникова https://zanauku.mipt.ru/2023/01/19/ot-nyutona-do-kovaleva/ (статья популярная, но там есть ссылки и на статьи с подробностями)

2.
на правой картинке — показано, к какому из комплексных корней сходится метод Ньютона, стартующий в данной точке

получается пример невозможной, казалось бы, ситуации: такие 4 области плоскости, что любая точка границы одной из них лежит на границе всех 4 («в любой точке границы граничат сразу 4 страны»)

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

мне кажется, впервые узнал про такое (как и про многое другое в математике) из рассказов Вити Клепцына — спасибо ему

***

ну или можно никаких объяснений не читать, просто посмотреть на красивые картинки… еще больше картинок у 3Blue1Brown: https://youtu.be/-RdOwhmqP5s

Посмотрим на то, для каких простых P уравнение X²=D разрешимо в остатках по модулю P.

На компьютерных экспериментах хорошо видно, что возникает удивительная периодичность — вот, например, для D=7 (для непростых модулей вместо ответа печатается точка):

x²=7 mod p
apparent (weak) period is 28
..++.-.+...-.-...-.+...-....
.+.+.....+...-.-...+.....+..
...+.-.....-...-.-.....-...+
.....-.......-...-.+...-.+..
.+.............-...+.....+.+
.........+.-.....-.....-...+
.....-.....-.-.........-.+..
.+.+...........-...........+
...+.-...+.....-.-.........+
.....-.....-.....-.+.....+..

Ответ всегда зависит только от остатка P mod 4D — это квадратичный закон взаимности (в форме, открытой Эйлером).

Можно было бы надеяться, что похожая история будет и в кубическом случае… но нет, ничего похожего. Простые P, для которых уравнение X³=2 (например) разрешимо mod P, выглядят хаотично, никакого периода нет. Единственная относительно видная закономерность — что доля простых, для которых X²=D разрешимо, всегда (ну если D не куб ) примерно 2/3.

x³=2 mod p
no apparent period
..++.+.-...+.-...+.-...+.....+.+.....-...+.+...+.....+.....+.-.....-...+.-.....-...+.....+.......-...+.-...+.+...+.............+...+.....+.-.........+.-.....+.....-...+.....+.....+.-.........+.-...+.-...........-...........+...+.+...+.....+.-.........+.....+.....+.....+.-.....+...+.+.........+.............(…)
+: 112 -: 56 rat: 0.67

Но можно посмотреть и на другие кубические многочлены — и иногда ответ снова определяется остатком P по какому-то модулю. Например:

x³-9x-9=0 mod p
apparent (weak) period is 9
..-+.-.-.
..-.-...+
.+...-...
..-.-....
.+...-.-.
..-.....+
.....-.-.
....-...+
.+.....-.
..-.....+
.......-.
.
+: 8 -: 17 rat: 0.32

x³-x²-4x-1=0 mod p
apparent (weak) period is 13
..--.+.-...-.
+...-.-...-..
...-.+.....-.
..-.-...+....
.+.....-.-...
..-...-.+....
.+...+.....-.
......-...-.+
...-.+...-...
..........-..
.+.....-.-...
......-.+....
.+.....-...-.
....-.....-.+
.........-.-.
..-.-
+: 14 -: 32 rat: 0.3

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

===

Это всё можно и в экселе делать — но в комментарии положу небольшую программу на питоне.

А для математического контекста — пусть здесь будет популярная лекция Фрэнка Калегари на ICM-2022: https://youtu.be/EDsK-8SBx-g (там как раз с таких вопросов все и начинается)

L-системы, или Как нарисовать снежинку и дерево

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

биолог Аристид Линденмайер в конце 1960-х придумал L-системы именно как модель для роста растений — и невозможно не написать про это два слова

сразу скажу, что есть отличная книжка Prusinkiewicz, Lindenmayer. The Algorithmic Beauty of Plants, http://algorithmicbotany.org/papers/#abop (спасибо Никите С. и Роме Г., советовавшим эту книжку) — и там все подробно step by step объясняется, лучшее мб вместо поста ее и читать )

===

3.
в предыдущем посте мы генерировали строки и смотрели на них как на описания путей (символ F означает «двигаться вперед», повороты кодируются символами + и -)

если хочется рисовать деревья, то полезно добавить какую-то возможность ветвления

для этого добавим такое правило для обработки символов [ и ]: после рисования фрагмента в квадратных скобках «черепашка» возвращается туда, где была в его начале

технически для этого в функцию turtleprint добавим

stack = []
…
if command == "[":
    stack.append((x,y,heading))
if command == "]":
    x,y,heading = stack.pop()

тогда дерево можно нарисовать в духе

seed="F"
rules={"F": "F[+F][-F]F"}

(стартуем с отрезка, на каждом шаге заменяем отрезок на цепочку из двух с растущими из середины в обе стороны ветками; на рисунке показано всего 2 шага)

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

seed="X"
rules={"F": "FF", "X": "F[+X][-X]FX"}

уже не так плохо, но слишком симметрично… перекосить можно так, скажем:

seed="X"
rules={"F": "FF", "X": "F-[[X]+X]+F[+FX]-X"}

(код целиком положу в комментарии)

***

(не слишком сложными) L-системами можно пользоваться и прямо в техе (в tikz есть для этого бибиотека):

\usetikzlibrary{lindenmayersystems}
\begin{tikzpicture}[rotate=65]
\draw[green!60!black] l-system[l-system={axiom=F,rule set={F -> F[+F]F[-F]}, order=4, angle=25,step=3pt}];
\end{tikzpicture}

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

===

вижу иногда, как люди постят generative art, где код умещается примерно в один твит, а на картинках облака летят, деревья растут, вселенные встают и рушатся, и всё такое

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

L-системы, или Как нарисовать снежинку и дерево

1.
начнем со слова A и дальше будем итерировать такие замены: A→AB, B→A (на каждой итерации одновременно заменяются все буквы по указанным правилам)

будут получаться всё большие отрезки интересного слова Фибоначчи — про него можно прочитать в статье Вити Клепцына в «Квантике» (№№5-6 за 2020 год)

такого рода правила (“L-системы”) легко реализовать — например, так:

def lsystem(seed,rules,iterations):
    for _ in range(iterations):
        seed = "".join(
            rules[letter] if letter in rules else letter
            for letter in seed)
    return seed

print(lsystem("A",{"A":"AB","B":"A"},5))

2.
такого рода правила порождают слова с фрактальной структурой… так что можно попробовать превратить их во фрактальные картинки

можно начать, например, со снежинки Коха: мы стартуем с правильного треугольника, а дальше на каждой итерации заменяем каждый отрезок на зубец («_ → _/\_»)

это можно закодировать так:

snowflake = lsystem("F++F++F",{"F":"F-F++F-F"},4)

где “F” означает «идем вперед и рисуем отрезок», а “+” и “-” — изменения направления

(про такое, кстати, тоже писали в «Квантике»: Валя Кириченко и Владлен Тиморин объясняли, как это связано со словом Туэ-Морса и проч.)

осталось это реально нарисовать — например, при помощи matplotlib:

import math
import matplotlib.pyplot as plt

def turtleprint(string,alpha0=0,alpha=90,step=1):
    fig = plt.figure()
    fig.set_size_inches(6,6)
    ax = fig.gca()
    ax.set_aspect('equal')

    x,y = 0,0
    heading = alpha0
    for command in string:
        if command == "F":
            x0,y0 = x,y
            x = x0 + step*math.cos(math.radians(heading))
            y = y0 + step*math.sin(math.radians(heading))
            plt.plot([x0,x],[y0,y],color='k')
        if command == "+":
            heading += alpha
        if command == "-":
            heading -= alpha

    plt.axis('off')
    plt.show()

для 4 итераций получается буквально картинка сверху

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

turtleprint(lsystem("-L",{"L":"LF+RFR+FL-F-LFLFL-FRFR+","R":"-LFLF+RFRFR+F+RF-LFL-FR"},4))

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

суммы косинусов и кубические уравнения

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

1.
напомню для начала квадратичный случай

пусть P простое число… и пускай вида 4k+1

квадратичная сумма Гаусса — это способ выразить √P через корни степени P из единицы с разными знаками (если показатель степени — квадратичный вычет mod P, то берем знак +1, для невычета — -1)

можно записать его в виде ∑cos(n² 2π/P)=√P

(например, для p=5: cos(2π/5)=(√5-1)/4, а вся сумма состоит из 4 равных косинусов и единицы)

2.
пускай теперь P=3k+1, а нас интересует тригонометрическая сумма G=∑cos(n³ 2π/P)

Галуа говорит нам, что G — корень кубического уравнения с рациональными коэффициентами, но какого конкретно?

полезно начать не с суммы G, а с кубической суммы Гаусса g

определяется она аналогично квадратичному случаю, но знак +1 теперь получают степени, являющиеся кубическими вычетами, еще два класса — знаки ω и ω² (ω — комплексный кубический корень из единицы)

(другими словами, g — линейная комбинация трех корней искомого уравнения, на которой группа Галуа действует просто умножением на ω)

снова |g|=√P, но теперь что-то хорошее получается при возведении не в квадрат, а в куб: g³=PП, где ПП' — разложение p на простые в Z[ω], а штрих обозначает комплексное сопряжение

3.
вернемся к тригонометрической сумме

G=g+g', поэтому (пользуясь описанными выше свойствами g)

G³=g³+g'³+3gg'(g+g')=P(П+П')+3PG,

т.е. G корень кубического уравнения x³-3Px-AP, где A=2ReП

чтобы получить уравнение буквально на суммы, которые рассматривались в предыдущем посте, нужно еще сделать замену x=-3t+1

4.
если хочется вычислить A без использования арифметики в Z[ω], то можно воспользоваться тем, что если P=ПП' в Z[ω], A=2ReП, то

4P=A²+27B²

такое представление единственно… с точностью до знаков; “наше” A — это то, которое дает остаткок 1 mod 3 (это вопрос про знак гауссовой суммы, который выше был заметен под ковер… там пусть и остается)

можно описать A и по-другому: количество решений по модулю P уравнения X³+Y³=1 равно P-2+A

5.
дискриминант уравнения x³-3Px-AP=0 как раз равен (27PB)²

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

уравнение x³=P, например, не подходит, а вот наше — вполне себе

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

6.
если P=9m²-3m+1, то 4P=(3m-2)²+27m². В этом случае свободный член уравнения на триг. суммы является полным кубом (числа m, собственно) и формула Рамануджана для суммы кубических корней из корней уравнения дает особенно простой ответ:

2S³ = 3³√(mP)-6m+1

например, для P=73 имеем m=3 и ответ ³√((3³√219-17)/2); для P=31 имеем m=-2 и опечатку в заметке в Мат. просвещении

на картинке сверху — тождества¹ из заметки С.Маркелова в Мат. просвещении, и там предлагается придумать обобщения

¹ там только есть опечатка… найдите

программа в комментариях — говорит, суммы каких косинусов надо взять для произвольного 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

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

различаются ли вычисленные приближенно рациональные и иррациональные числа?

про математику здесь хорошо написано в самом начале книги «Математический дивертисмент» Табачникова и Фукса (которую вообще всячески рекомендую)

а я только покажу программу, которая различает рациональные и иррациональные числа: по числу, например, 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)

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

что делать с ответами? (продолжение)

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

что это конкретно за прямая? можно не пытаться подобрать коэффициенты на глаз и т.п., а воспользоваться функцией экселя LINEST

если log P(2N) ≈ kN+b, то P(2N) ≈ c⋅λ^N (для λ=10^N, c=10^b)

хорошее ли получается приближение? вот фиолетовым цветом внизу таблицы показаны эти самые приближения c⋅λ^N, а рыжим цветом в последней строчке — величина ошибки

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

эксель в комментариях

===

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

что делать с ответами?

начнем давайте с модельной задачи… вот хотя бы с теми же доминошками: посчитаем количество P(N) разбиений на доминошки прямоугольника 3×N

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

хорошо, сделали файл, в котором посчитаны сколько-то первых значений P(N) — и что дальше?

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

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

когда целые числа быстро растут, то совсем грубо можно смотреть просто на количество цифр в P(N) — ну или если не на таком детском языке, то можно говорить про логарифмическую шкалу для оси y

если это сделать и прострить график не P(2N), а log P(2N), то… ну это практически как маленькое чудо выглядит (и на школьников, вроде, производит некоторое впечатление)

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

Дима Швецов напомнил про такой сюжет

Квантик пытается экспериментально изучить сумму обратных простых

Ноутик воспользовался освоенной им недавно библиотекой 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

дайджест комментариев: разбиения на доминошки

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

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))

как делить числа?

пусть у нас есть числа 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/

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

суммы двух квадратов обсуждают часто, а родственные утверждения — реже

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

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

https://sites.math.rutgers.edu/~zeilberg/expmath/

есть и видеозаписи

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

² см. https://t.me/mathtabletalks/733 и далее

===

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

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

можно попробовать начать с какого-то фиксированного разбиения (например, со всеми горизонтальными доминошками), а дальше применять случайные флипы: выбираем случайный квадрат 2×2 в фигуре, если он разбит на две доминошки, делаем флип (= <-> ||)

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

на первой картинке для квадрата 20×20 без угла сделано 100 флипов и видно, что мы пока от начального разбиения недалеко ушли; на второй картинке 10000 флиплов — и это выглядит уже довольно хаотично

вроде все ожидаемо, ничего особо интересного

...

попробуем теперь то же самое сделать с ацтекским брильянтом порядка 20… ну первый сюрприз, что после 10 тыс. флипов не происходит почти ничего… это, впрочем, не удивительно, в начале даже подходящих квадратов совсем мало…

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

===

способ хранения разбиения из предыдущего поста оказался неудобным тем, что значок не говорит, в какую сторону продолжается покрывающая доминошку клетка — стал хранить в каждой клетке указатель на парную к ней (не в программистском смысле указатель, а один из символов →←↓↑ )

картинки рисую matplotlib'ом (код в комментариях)

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