…в частности, сгенерировал случайное замощение брильянта порядка 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

1.
понравился прием с reshape — покажу на примере: пускай мы интересуемся, скажем, представлением чисел в виде суммы квадратов

import numpy as np

def sum_squares(N,k):
    counts1 = np.zeros(N*N+1,dtype=int)
    counts1[0] = 1
    np.add.at(counts1,[n*n for n in range(1,N+1)],2)
    counts = counts1.copy()
    for _ in range(k-1):
        counts = np.convolve(counts,counts1)
    return counts[:N*N+1]

(такой подсчет сумм сверткой — это примерно как возводить производящую функцию квадратов в степени, только без “развлечений” с символическими вычислениями)

тогда можно печатать ответы по сколько-то в строку в духе

print(np.reshape(sum_squares(20,3)[:240],(-1,8)))

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


2.
для сумм 4, скажем, квадратов ничего из reshape вроде бы не видно… зато помогает посмотреть на значения для простых чисел — на бумажке это немного утомительно, а компьютерно легко (все числа там получаются кратными 8, так что сразу на 8 разделим):

counts = sum_squares(N:=10,4)
for p in np.nonzero(is_prime(N*N))[0]:
    print(f"{p}: {counts[p]//8}")

какая возникает гипотеза? )

а что будет для p², p³? (спойлер: ответ равен количеству одномерных подпростравнств в 2-, 3-, 4-мерных пространствах над Z/p — это приятно понимать при помощи кватернионов)

(upd)
заодно посчитал количество представлений простых в виде сумму 8 квадратов — и тут легко угадал и неизвестный мне заранее ответ: 16(p³+1) — мб надо подумать, как это доказать

бонусная картинка в продолжение предыдущего: если простое число представимо в виде суммы а) 2; б) 3; в) 4; г) 8 квадратов, то сколько таких представлений

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

(там еще добавлены для мастштаба две линии типа с√n — но на самом деле максимумы растут примерно как √n log log n, насколько понимаю)

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

Г.Б.Шабат попросил продемонстрировать это на примере уравнения x^3-8x^2+5x+1=0 (дискриминант 7^4) — оказалось, как нередко бывает, что между «в принципе понимаю как делать» и «могу сделать на практике» есть некоторый зазор (но постепеннно все получилось)

1.
легко попросить sympy решить кубическое уравнение:

from sympy import *

x = Symbol('x')
P = x**3 - 8*x**2 + 5*x + 1
theta0, theta1, theta2 = solve(P,x)

но после print(theta0) мы увидим длинную формулу с двумя кубическими корнями из комплексных чисел, прямо из нее непонятно как извлечь ответ с косинусами

2.
чтобы решать кубическое уравнение, удобно сделать преобразование Фурье и смотреть не на корни, а на такие их линейные комбинации («резольвенты Лагранжа»):

omega = exp(2*pi*I/3)
lambda0 = theta0 + theta1 + theta2
lambda1 = theta0 + omega*theta1 + omega**2*theta2
lambda2 = theta0 + omega**2*theta1 + omega*theta2

чем это полезно? тем, что когда группа Галуа переставляет тэты по циклу — lambda1 и lambda2 домножаются на степени ω, соотвественно их кубы неподвижны, т.е. лежат в Q[ω]

собственно, это в формуле для корней кубического уравнения мы и видим: theta0=(lambda0+lambda1+lambda2)/3, и первое слагаемое — коэффициент нашего уравнения (в данном случае, 8), а вторые два — кубические корни каких-то крокодилов, лежащих в Q[ω]

3.
как найти нужного крокодила? print(lambda1**3) дает невразумительное выражение на несколько строк, а хотелось бы видеть конкретный ответ a+bω с рациональными (на самом деле, даже целыми) a и b… спросим sympy получше:

print(minimal_polynomial(lambda1,x,domain=QQ))

x**6 - 637*x**3 + 117649

число 117649 может пугать (если не сразу заметили, что это 7^6) — но в любом случае, мы получили просто квадратное уравнение (на куб lambda1), его можно решить и получить -7²(1+3ω)²

(можно на всякий случай проверить: print(N(lambda1**3+7**2*(1+3*omega)**2)) дает более-менее ноль)

если читали https://www.group-telegram.com/compmathweekly.com/23 — можно подумать, как теперь получить конкретный ответ )

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

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

4.
итак, корни кубического уравнения выражаются через лямбды, и в нашем случае λ₁³=-7²(1+3ω)²

теорема Кронекера-Вебера (whatever it means) гарантирует, что λ₁ можно выразить через корни из единицы, но как конкретно это сделать?

выше обсуждалось, что если π=a+bω простое с нормой p, то кубическая сумма Гаусса по модулю p — это примерно кубический корень из pπ

спускаясь от общих разговоров к конкретным вычислениям: N(1+3ω)=7, мы хотели бы извлечь кубический корень из 7(1+3ω), попробуем взять эту самую сумму Гаусса:

zeta = exp(2*pi*I/7)
g = sum(omega**k * zeta**(3**k) for k in range(6))
print("-g^6 vs lambda1^3:",N(-g**6),N(lambda1**3))
print(“-g^2 vs lambda1:”,N(-g**2),N(lambda1))

(число 3 в определении g — это первообразный корень mod 7; отметим, что g — это, как и было обещано, линейная комбинация каких-то корней из единицы)

запустив программу можно увидеть, что g³=7(1+3ω) (могли бы получить сопряженное и/или какой-то знак), но -g² — не тот кубический корень, который нам нужен; нетрудно перебрать варианты, правильный ответ λ₁=-ωg²

l1 = -omega*g**2
print("l1 vs lambda1:",N(l1),N(lambda1))

5.
осталось превратить это в выражение для корней

по определению лямбд, θ₀=(λ₀+λ₁+λ₂)/3; в нашем случае λ₀=8, а λ₁ и λ₂ комплексно сопряжены, так что можно написать

t0 = (2*re(expand(l1))+8)/3
print("t0 vs theta0",N(t0),N(theta0))
print(t0)

так как всё уже выражено через корни из единицы, после взятия вещественной части sympy сразу даст ответ, записанный через косинусы:

-10*cos(2*pi/21)/3 - 4*cos(2*pi/7) - 2*cos(4*pi/21) - 4*cos(8*pi/21)/3 - 2*cos(10*pi/21) + 4*cos(pi/21)/3 + 8*cos(pi/7)/3 + 10*cos(5*pi/21)/3 + 8/3

6.
можно выразить и остальные корни

t0 = (2*re(expand(l1))+8)/3
t1 = (2*re(expand(l1*omega**2))+8)/3
t2 = (2*re(expand(l1*omega))+8)/3

как связаны получающиеся выражения?

вот они (слегка переписанные — в одинаковом порядке и т.п.):

t0 = 8/3
- (10/3)*(cos(2*pi/21) + cos(16*pi/21))
- (8/3)*cos(6*pi/7) - 4*cos(2*pi/7)
- 2*(cos(4*pi/21) + cos(10*pi/21))
- (4/3)*(cos(8*pi/21) + cos(20*pi/21))

t1 = 8/3
- (10/3)*(cos(4*pi/21) + cos(10*pi/21))
- (8/3)*cos(2*pi/7) - 4*cos(4*pi/7)
- 2*(cos(8*pi/21) + cos(20*pi/21))
- (4/3)*(cos(16*pi/21) + cos(2*pi/21))

t2 = 8/3
- (10/3)*(cos(8*pi/21) + cos(20*pi/21))
- (8/3)*cos(4*pi/7) - 4*cos(6*pi/7)
- 2*(cos(16*pi/21) + cos(2*pi/21))
- (4/3)*(cos(10*pi/21) + cos(4*pi/21))

— видно, что следующее выражение можно получить, увеличив в предыдущем все аргументы вдвое

нетрудно и без выписывания этих сумм понять, почему так должно быть: автоморфизм ζ→ζ², ω→ω² поля Q(ζ,ω) домножает g на кубический характер двойки mod 7, т.е. на ω², т.е. λ₁=-ωg²→-ω²(ω²g)²=λ₁/ω etc — то есть действительно сдвигает корни на один по циклу

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

наверное от группы можно было и плясать — перебирать элементы порядка 3 в Z/2×Z/6, для каждого варианта писать свои суммы, считать соотв. мин. многочлены, то-сё… не продумывал

для недавно присоединившихся — некоторые старые посты:

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

Все знают, как найти сумму 1+2+…+N. Некоторые помнят формулу для суммы 1²+2²+…+N². А как найти формулу для суммы 4-х, например, степеней?

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

Сделаем табличку для первых нескольких степеней первых нескольких чисел, а также для первых сумм 1^5+…+N^5. Мы надеемся, что для этой суммы есть полиномиальная формула — то есть что колонка сумм линейно выражается через остальные. Чтобы искать (приближенно) линейную зависимость, можно воспользоваться уже знакомой по https://www.group-telegram.com/compmathweekly.com/20 функцией LINEST.

Получаем (гипотетический) ответ
1^4 + … + N^4 = N^5/5 + N^4/2 + N^3/3 - N/30

Можно добавить еще одну колонку и проверить, что для небольших N эта формула дает правильные ответы. Можно взять бумажку и ручку и доказать всё по индукции.

А можно подумать еще немножко и понять, что так найденный ответ автоматически будет правильным… — писал про это немного в Кванте, см. теорему 1 в https://dev.mccme.ru/~merzon/pscache/bernoulli-howto-pre.pdf (а в конце, кстати, там есть и код на питоне).

коллега Шатохин показал забавное развлечение (в т.ч. для не супер продвинутых участников):

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

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

(ну например: синяя диагональ возникает потому, что N² = 0 mod N и вообще (N+k)² mod N маленькое число, когда k маленькое… и т.д. и т.п.)

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

начинал уже ( https://www.group-telegram.com/compmathweekly.com/45 ) разговор про подсчет количеств решений mod p

базовый пример здесь — плоские кривые: пишем уравнение на x и y с целыми коэффициентами и смотрим как растет количество решений mod p с ростом p

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

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

дальше последовательность выглядит как «один, два, много» — кубические кривые уже скрывают бесконечную сложность… но чтобы с чего-то начать:

если мы смотрим на число N(p) решений y²=x³+ax²+bx+c mod p (и всё гладко, что бы это ни значило… напр., кривая y²=x³ не подходит), то можно ожидать, что правая часть примерно с одинаковой вероятностью квадратичный вычет и квадратичный невычет… и если воспринимать здесь идею про случайность всерьез, то можно ожидать, что |N(p)-p| имеет порядок примерно √p

из (доказанных) гипотез Вейля следует, что |N(p)-p|⩽2√p, в частности, N(p)/p→1… а дальше можно посмотреть на произведение N(p)/p (по p⩽x) — и ожидается, что эта штука растет примерно как log(x)^r, где r — ранг нашей кривой (рациональные точки на кривой образуют коммутативную группу, речь идет про ее ранг)

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

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

#видеозаписи #ТяжелыйПонедельник

Открываем новую видеозапись выступления:
Иван Дианов — Алгоритмическая графика: JS для дизайна и искусства

YouTube | VK Видео

1.
сколько коник проходят через 5 фиксированных точек общего положения на плоскости?

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

в данном случае, коника в P^2 задается 6 коэффициентами ее уравнения, все коники образуют проективное пространство P^5

каждое условие “проходить через точку” — гиперплоскость в этом пространстве

5 гиперплоскостей общего положения в 5-мерном пространстве пересекаются ровно по одной точке

(упражнение к части 1: сколько рациональных кубик проходит через 8 точек общего положения на плоскости?)


2.
продолжим разминку: сколько прямых пересекают 4 фиксированные прямые общего положения в пространстве?

(дальше будут когомологии и характеристические классы, простите… если это пугает, посмотрите вместо этого на Шуховские башни — https://book.etudes.ru/articles/shuhov/ — и поймите ответ на вопрос про прямые без всякий когомологий)

теперь нас интересует пространство PGr(1,3) всевозможных прямых в P^3 (оно же Gr(2,4), пространство 2-мерных линейных подпространств в 4-мерном); каждое условие «пересекаться с данной прямой» задает в этом пространстве гиперповерхность — и мы хотим посчитать число точек пересечения 4 таких гиперповерхностей

в проективном пространстве теорема Безу говорит, что (в хорошей ситуации) достаточно перемножить степени интересующих нас гиперповерхностей

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

а классы, которые мы будем там перемножать, — это обычно хар. классы каких-то естественных расслоений, в данном случае — все прямые, пересекающие данную, представляют c_1(E^*), где E тавтологическое расслоение (упражнение для тех, кому понятна формулировка: убедить себя в этом)

в качестве (аддитивного) базиса в H(Gr(n,n+m)) можно взять многочлены Шура, которые нумеруются диаграммами Юнга внутри прямоугольника n×m



то есть в принципе можно пропустить все разговоры про когомологии и т.п. как страшный сон — операционально всё сводится к элементарной алгебре многочленов и/или элементарной комбинаторике диаграмм Юнга

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

3.
итак, симметрические функции в sage-101:

Sym = SymmetricFunctions(QQ)
Sym.inject_shorthands()

теперь можно сказать, например,

s[1].expand(2)

и посмотреть как выглядит многочлен Шура для одной клеточки (как сумма всех переменных — в данном случае, двух: x0+x1)

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

R = PolynomialRing(QQ,'x',2)
R.inject_variables()
s.from_polynomial(x0^2+x1^2)

-s[1, 1] + s[2]

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

s[1]^4

s[1, 1, 1, 1] + 3*s[2, 1, 1] + 2*s[2, 2] + 3*s[3, 1] + s[4]

(ну конкретно это вычисление можно сделать и в уме — если помнить, что умножение на s[1] всевозможными способами добавляет клетку… но не будем отвлекаться)


4.
когомологии Gr(n,n+m) — это фактор кольца симметрических функций по многочленам Шура диаграмм, не помещающихся в прямоугольник n×m

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

базис в H(Gr(n,n+m))— это многочлены Шура от диаграмм, помещающихся в прямоугольник m×n, причем многочлен Шура столбца длины k — это k-й класс Черна двойственного тавтологическому расслоения E* (т.е. виртуальные переменные, от которых все многочлены симметрические, можно отождествлять с корнями Черна для E*), а единственная образующая максимальной степени, многочлен Шура для прямоугольника n×m — это класс точки в старших когомологиях

то есть, например, в когомологиях PGr(1,3)=Gr(2,4) мы посчитали, что c_1(E*)^4 = 2[pt]

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

и на разные другие вычисления с многочленами Шура можно смотреть геометрически в духе Шуберта… скажем, s[1]²=s[1,1]+s[2] — это утверждение о том, что гомологически [прямые, пересекающие две фиксированные]=[прямые в фиксированной плоскости]+[прямые через фиксированную точку] — и геометрически это можно увидеть, если две фиксированные прямые сделать пересекающимися


5.
попробуем что-то немного посерьезнее: посчитаем количество прямых на кубической поверхности в P^3 общего положения

уравнение кубической поверхности дает сечение расслоения Sym^3 E* на том же PGr(1,3); нули этого сечения соответствуют прямым, ограничение на которые нашего уравнения тождественно равно 0, — т.е. в точности прямые, лежащие в нашей поверхности

итак, прямые на (типичной) кубической поверхности представляют класс Эйлера (он же старший класс Черна) расслоения Sym^3 E*

если корни Черна расслоения E* — x0 и x1, то у расслоения Sym^3 E* это 3x0, 2x0+x1, x0+2x1, 3x1; осталось их перемножить и понять, кто получился в терминах образующих нашего кольца — то есть операционально достаточно спросить

s.from_polynomial( 3*x0*(2*x0+x1)*(x0+2*x1)*3*x1 )

и в ответе

81*s[1, 1, 1, 1] - 63*s[2, 1, 1] + 27*s[2, 2] + 18*s[3, 1]

посмотреть на коэффициент при s[2,2] — вот они, 27 прямых

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

аналогичным образом можно посчитать количество прямых на типичной гиперповерхности степени 5 в P^4

s.from_polynomial( 5*x0*(4*x0+x1)*(3*x0+2*x1)*(2*x0+3*x1)*(x0+4*x1)*5*x1 )

15625*s[1, 1, 1, 1, 1, 1] … + 2875*s[3, 3] - …

(но конечно здесь [снова] заметено под ковер, почему мы считаем то, что надо: пересечения трансверсальны и т.п.)

давайте на этот раз что-нибудь более приземленное

знаете ли вы, чему равна сумма

a²/(a-b)(a-c) + b²/(b-c)(b-a) + c²/(c-a)(c-b)?

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

предлагаю поставить реакцию
«ok», если знаете,
«✓», если сделали упражнение


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

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

def sum1(n,k):
    R = PolynomialRing(QQ, 'x', n+1)
    x = R.gens()
    return sum(
            x[i]^k / prod(x[i] - x[j] for j in range(n+1) if j != i)
            for i in range(n+1)
        )

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


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

P.S.
на самом деле, этот сюжет связан из предыдущим: например, самое верхнее выражение — это подсчет прямых через две точки плоскости при помощи формулы локализации Атьи-Ботта на P²

5 лет назад был отличный математический флэшмоб #12equations — в т.ч. как раз в этот день пять лет назад писал про тангенс

…Очень люблю цитату из интервью Гельфанда, «я считал, что есть две математики — алгебраическая и геометрическая, и что геометрическая математика принципиально “трансцендентна” для алгебраической. (…) Когда я обнаружил, что синус можно записать алгебраически в виде ряда, барьер обрушился, математика стала единой».

В отличие от синуса и косинуса ряд для тангенса на первый взгляд выглядит хаотично:
x+x^3/3+2x^5/15+17x^7/315+…
Оказывается, за этим хаосом скрывается отличная комбинаторика…

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

x = var('x')
taylor = tan(x).taylor(x,0,21)
E = [ taylor.coefficient(x,n)*factorial(n) for n in range(1,22,2) ]
print(E)

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

или вот такой листок для курса Е.Смирнова про этот сюжет делал: https://dev.mccme.ru/~merzon/ium-combi/combi20-05-bernoulli.pdf

в продолжение, если угодно, темы про L-системы ( https://www.group-telegram.com/compmathweekly.com/25 & https://www.group-telegram.com/compmathweekly.com/33 ) —

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

(в качестве интермедии — никакой математики, просто картинки)

геометрия пи в экселе

выше обсуждались разные способы вычисления π, в т.ч. 1-1/3+1/5-1/7+… = π/4 ( https://www.group-telegram.com/compmathweekly.com/42 )

а можно ли как-то несложно поставить компьютерный эксперимент, где π возникает сразу «геометрически», без анализа?

π — это площадь круга радиуса 1, поэтому среди точек квадрата 0⩽X,Y⩽1 доля точек, попадающих в (четвертинку) круга радиуса 1 с центром в начале координат, равна π/4

теперь можно вычислять π геометрически хоть в экселе — по методу Монте-Карло: берем 100500 строк, в каждой пишем по паре случайных чисел X и Y, считаем в скольких строках X²+Y²⩽1, полученную долю умножаем на 4

конечно, способ не особо эффективный — для 1000 строк получается приближение уровня π≈3,1

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

P.S.
сумма из начала поста на самом деле получается, если думать про те же площади — но брать не случайные точки, а все точки достаточно мелкой сетки… про это есть отличный ролик 3Blue1Brown, если кто не видел, https://youtu.be/NaL_Cb42WyY

6.
попробуем еще чуть-чуть продвинуться в подсчетах кривых

в прошлый раз ( https://www.group-telegram.com/compmathweekly.com/72 ) посчитали количество прямых на гиперповерхности степени 5 в P^4

теперь попробуем посчитать на ней же коники

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

каждая коника в P^4 лежит в какой-то плоскости и высекается там квадратным уравнением... то есть пространство C всех коник — это P(Sym^2 E), где E — двойственное¹ тавтологическому трехмерное расслоение на Gr(3,5)=PGr(2,4)

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

уравнение f нашей квинтики задает сечение расслоения Sym^5 E, а условие «данная коника лежит на квинтике» значит, что f [в ограничении на плоскость коники] = (уравнение коники)⋅(что-то кубическое) — т.е. нас интересуют точки, для которых сечение попадает в подрасслоение (Sym^3 E)⊗O(-1) — другими словами, нас интересует старший класс Черна фактора F := Sym^5 E / (Sym^3 E)(-1)

(проверка: dim C=11, т.к. оно расслаивается над 6-мерным грассманианом с 5-мерным слоем; dim F=7⋅6/2-5⋅4/2=11 — размерности сходятся)

осталось реализовать компьютерное вычисление интеграла c_11(F) по C

(
прямолинейно-алгебраическое вычисление вот какое:

пусть x0, x1, x2 — корни Черна E, тогда у Sym^n E корни Черна — это i⋅x0+j⋅x1+k⋅x2 | i+j+k=n

H(C) = H(Gr)[t] / prod(t-x), где t — класс Черна O(1), а произведение берется по корням Ч. для Sym^2 E

нам нужно взять произведение 1+x по корням Ч. для Sym^5 E, разделить на произведение 1+x-t по корням Ч. для Sym^3 E, разложить t в степени 6 и выше по соотношению выше, найти коэффициент при s[2,2,2]⋅t^5

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