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

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


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/

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


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

справа можно видеть фрагмент квазипериодического замощения плоскости

в нем участвуют равнобедренные треугольники с углами при вершине π/5 (красные, «A») и 3π/5 (синие, «B»)

они замечательны тем, что A можно разбить на уменьшенные копии A,B,A, ну а B можно разбить на уменьшенные копии A,B — и если начать с А и итерировать такие замены, то можно думать, что мы собираем из треугольников A и B всё большую копию треугольника¹ A (в левой половинке картинки — первая пара итераций)

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

а вечером подумал, что это просто L-система — только параметрическая: кроме буквы A/B нужно помнить, как именно треугольник расположен на плоскости (и правила замены эти параметры должны правильно менять) — так что можно быстренько реализовать

¹ а чтобы получить замощение плоскости, можно, скажем, стартовать с 10 треугольников A с общей вершиной



положение треугольника решил хранить в виде пары комплексных чисел² (преобразования z→az+b, переводящего эталонный треугольник в наш) и написал такой шаг для получающейся параметрической l-системы:

phi = (math.sqrt(5)+1)/2
rot = math.cos(math.pi/5)+math.sin(math.pi/5)*1j

def step(state):
    for atom in state:
        c, a, b = atom
        if c=='A':
            yield ('A',a,b*phi)
            yield ('B',a*(rot**4),(a+b)*phi)
            yield ('A',a*(rot**3),(a+b)*phi)
        if c=='B':
            yield ('A',a,b*phi)
            yield ('B',a*(rot**4),(a+b)*phi)

state = [('A',rot**i,0+0j) for i in range(10)]
for _ in range(6):
    state = step(state)

по сути на этом все! — остается только дописать код для рисования треугольничков… ну программа целиком будет в комментариях

² уже засомневался, так ли это удачно — потому что для настоящей мозаики Пенроуза треугольники полезно и переворачивать

This may depend on your specific area, but in bijective combinatorics you have to believe in miracles! Otherwise you can never fully appreciate prior work, and can never let yourself loose enough to discover new miracles. To give just one example, the RSK correspondence is definitely on everyone’s the top ten list of miracles in the area. By now there are at least half a dozen ways to understand and explain it, but I still consider RSK to be a miracle.

(из недавнего поста Игоря Пака)

про RSK — невозможно, кажется, поспорить

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

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

в т.ч. длина первой строки P — как у самой длинной (нестрого) возрастающей подпоследовательности, а длина первого столбца — как у самой длинной (строго) убывающая подпоследовательности (бесплатное следствие — теорема Эрдеша-Секереша: в последовательности длины nm+1 либо первая длиннее n, либо последняя длиннее m)

заодно можно компьютерно убедиться в том, что для случайной перестановки N элементов самая длинная возрастающая подпоследовательность имеет длину ~2√N (как доказали Вершик-Керов, Логан-Шепп)

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

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

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

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

попросил chatgpt сконвертировать реализацию RSK из прошлого поста в javanoscript, поправил минимально оформление — и пожалуйста, вот веб-версия:

https://dev.mccme.ru/~merzon/compmath/RSK.html

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

Ваня Яковлев напомнил задачу о разборчивой невесте — вот формулировка из брошюры С.М.Гусейна-Заде:

В некотором царстве, в некотором государстве пришло время принцессе выбирать себе жениха. В назначенный день явились 1000 царевичей. Их построили в очередь в случайном порядке и стали по одному приглашать к принцессе. Про любых двух претендентов принцесса, познакомившись с ними, может сказать, какой из них лучше. Познакомившись с претендентом, принцесса может либо принять предложение(и тогда выбор сделан навсегда), либо отвергнуть его (и тогда претендент потерян: царевичи гордые и не возвращаются).Какой стратегии должна придерживаться принцесса, чтобы с наибольшей вероятностью выбрать лучшего?

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

но остается вопрос, как выбрать момент s (долю от 0 до 1) перехода от первой стадии ко второй — и здесь вполне можно поставить компьютерный эксперимент:

def concrete_test(s,sigma):
    N = len(sigma)
    n = math.floor(N*s)
    m = max(sigma[:n]) if n>0 else -1
    # пропустили долю s кандидатов
    # и ждем первого, лучше всех предыдущих
    while n<N and sigma[n]<=m: n += 1
    return n<N and sigma[n]==N-1 # нашли максимум?

def random_tests(s,N,T):
    sigma = np.array(range(N))
    wins = 0
    for _ in range(T):
        random.shuffle(sigma)
        if concrete_test(s,sigma): wins += 1
    return wins/T

попросил и из этого chatgpt сделать html-версию и мгновенно… получил программу с неправильной генерацией случайной перестановки

но вроде всё удалось исправить — и по ссылке https://dev.mccme.ru/~merzon/compmath/bride.html можно выбрать момент s и посмотреть, насколько успешна ваша стратегия

для меня сюрпризом стало, что в реальности разница между теоретически лучшим вариантом (ну почитайте Гусейна-Заде!) и наивным s=50% (или, скажем, s=25%) ничтожна

(тут, кстати, легко грубо прикинуть: вероятность того, что лучший кандидат во второй половине, а second best в первой, — уже 1/4, так что s=1/2 приводит к успеху с вероятностью больше 25%)

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

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

их оказывается столько же, сколько матриц из целых неотрицательных чисел, если число в клетке (i, j) считать с весом i+j-1 (так что, например, производящая функция для всех трехмерных диаграмм равна Π(1-q^i)^{-i})

это может показаться смутно гармоничным (если в плоском разбиении на клетке (i,j) лежит кубик, то и на клетках (i',j) и (i,j') при всех i'<i и j'<j лежат кубики… и всего таких кубиков минимум i+j-1…), но если пытаться придумать конкретную биекцию, возникает ощущение магии

это тоже магия соответствия RSK


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

реализация — в комментариях


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

плюнул и за 10 минут написал matplotlib-рисование кубиков (результат на картинке)… ну может еще и сделаю веб-версию, не знаю

***

у нас с Мишей Берштейном есть текст в МатПросвещении, где рассказывается в том числе про подсчет трехмерных диаграмм Юнга — но там никакого RSK как раз нет, другое есть: https://www.mathnet.ru/rus/mp823

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

сколькими способами число N можно представить в виде суммы?

(разбиения 1+2 и 2+1 и т.п. считаются одинаковыми)

начать эксперименты можно и с примитивного перебора — например, такого:

def partitions(n,l=None):
    # все разбиения числа n
    # (на числа не превосходяшие l)
    if l is None: l=n
    if n==0: yield []
    if n>0:
        for i in range(l,0,-1):
            # рекурсивный перебор по максимальному элементу
            yield from map(lambda partition: [i]+partition,
                           partitions(n-i,i))

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

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

def count_partitions_upto(N):
    partitions = [[0 for _ in range(N+1)] for _ in range(N+1)]
    for l in range(N+1):
        partitions[0][l] = 1
    for n in range(1,N+1):
        for l in range(1,N+1):
            for i in range(1,l+1):
                partitions[n][l] += partitions[n-i][i]
    return partitions

получаются какие-то числа… что дальше с ними делать? как обсуждалось выше, полезно бывает посмотреть на график — и, так как p(N) быстро растет, лучше смотреть на график log p(N)

картинки (и еще пара слов) в комментариях

дальше мб будет что-то про p(N) по разным модулям

в дополнение к предыдущему — количества p(n) разбиений — по разным модулям

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

спойлер:
p(5n+4) делится на 5
p(7n+5) делится на 7
p(11n+6) делится на 11

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

в качестве мат. контекста можно прочитать короткую заметку Ramanujan’s congruences and Dyson’s crank (G.Andrews, K.Ono)

***

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

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

***

заодно решил посмотреть на количества разбиений с разными рангами (контекст e.g. в заметке выше) — вот код для этого:

def partitions_upto(N):
    partitions = [[] for _ in range(N+1)]
    partitions[0] = [(".",0,0)]
    for n in range(1,N+1):
        for i in range(n,0,-1):
            for pstr,pmax,plen in partitions[n-i]:
                if i>=pmax:
                    partitions[n].append((str(i) if pstr=="." else f"{str(i)}+{pstr}", i, plen+1))
    return partitions

mod = 7
k = 4
partitions = partitions_upto(mod*k)
print(f"mod = {mod}")
for n in range(mod*k+1):
    ans = [0 for _ in range(mod)]
    for pstr,pmax,plen in partitions[n]:
        ans[(pmax-plen) % mod] += 1
    print(f"{n} — {len(partitions[n])} ({len(partitions[n]) % mod}):", *ans)

о том, как (не) стоит считать π, и о магии им. Эйлера


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

если взять arctg(x) = x - x³/3 + x⁵/5 - …, подставить¹ x=1, то получается замечательно простая формула Лейбница π/4 = 1 - 1/3 + 1/5 - 1/7 + … — чего еще желать?..

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

в прошлый раз для сотни правильных знаков хватило 4 итераций, сейчас попробуем взять… ну, скажем, 1000 членов:

from mpmath import *

n = 1000
mypi = mpf(0)
s = mpf(4)
for k in range(n):
    mypi += s / mpf(2*k+1)
    s = -s

print(n,nstr(mypi,10),
      "diff:",nstr(pi-mypi,3,min_fixed=1))

1000 3.140592654 diff: 1.0e-3

видно, что погрешность ~1/1000… и вообще, можно понять, что для N слагаемых такой способ дает погрешность ~1/N — то есть чтобы получить 100 знаков нам нужно сложить ~10^(100) слагаемых (!!!)


2.
кажется, что это совсем тупиковая ветвь… и все же, не будем унывать сразу

посмотрим, скажем, на сумму 5 млн членов

погрешность получается ожидаемая, порядка 1/N, т.е. шесть цифр после запятой верны, а седьмая уже нет (получается 3,1415924… вместо 3,1415926…)

но то, что видно дальше, иначе как магией не назовешь: например, следующая дюжина цифр снова правильная!

from mpmath import *
mp.dps = 100

n = 2500000
mypi = mpf(0)
for k in range(1,n+1):
    mypi += mpf(8) / ((4*k-3)*(4*k-1))

for (s1,s2) in zip(nstr(mypi,60),nstr(pi,60)):
    print(s1 if s1==s2 else s1+'\u0332',end='')
print()

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

3.14159245358979323846464338327950278419716939938730582097494

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


продолжение ниже

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

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

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

просуммируем первые N членов… насколько мы далеки от правильного ответа? сумма оставшихся членов примерно равна интегралу 1/x³ от N до ∞, т.е. ~1/2N²

то есть даже для N=100500 мы получаем всего 10 правильных цифр, а до 100 правильных цифр дойти практически невозможно (требуется ~10^{50} слагаемых)…


…но подождите! «ваш прогноз выполняется с вероятностью 40%? так говорите все наоборот, и будет выполняться с вероятностью 60%» — раз мы знаем, насколько наш ответ далек от правильного, так давайте соответствующим образом скорректируем ответ!

и действительно, для N=100500 сумма первых N членов +1/2N² дает на 5 больше правильных цифр


можно ли оценить сумму (в данном случае, хвост ряда) точнее, чем просто интегралом? оценка «интеграл ≈ сумма» возникает при приближении площади под графиком площадью столбиков — а ясно, что лучше заменить прямоугольники трапециями

если вдуматься, это учит, что точность повысится, если вычесть половину N-го члена

на самом деле, +1/2N²-1/2N³ — это начало серии поправок, а возникающие коэффициенты — это (примерно) числа Бернулли… писал про этот сюжет (формулу Эйлера-Маклорена и числа Б.) в статье «Алгебра, геометрия и анализ сумм степеней последовательных чисел», https://www.mathnet.ru/rus/mp882

***

посмотрим, как это работает на практике

возьмем, скажем, 57 поправочных членов… а сумму не обязательно даже считать до N=100500, ограничимся N=500 — и уже этого хватает для 100 правильных цифр после запятой!

from mpmath import *
mp.dps = 102

from textwrap import wrap

def bernoulli_upto(m):
    ans = [mpf(1)]
    for k in range(1,m+1):
        b = mpf(0)
        for i in range(k):
            b -= ans[i]*binomial(k+1,i) / (k+1)
        ans.append(chop(b,tol=1e-10))
    return ans

def improve(pre,N,m):
    bernoulli = bernoulli_upto(m-1)
    for k in range(1,len(bernoulli)+1):
        pre += (bernoulli[k-1]*k) / (2*N**(k+1))
    return pre

def partial_sum(N):
    ans = mpf(0)
    for k in range(1,N+1):
        ans += mpf(1) / k**3
    return ans

def mprint(num,prec):
    intp, frcp = nstr(num,prec).split(".")
    print(f"{intp},",*wrap(frcp,width=5))

#N = 100500
#ans = partial_sum(N)
#mprint(ans,21) # всего 10 правильных цифр

N = 500
ans = partial_sum(N)
mprint(improve(ans,N,57),101) # все 100 цифр правильные!

***

так что, может вообще зафиксировать N, а только увеличивать количество поправочных членов (аргумент m функции improve)? степени N экспоненциально убывают, ну какие-то коэффициенты стоят при них…

проблема в том, что если начало последовательности Бернулли выглядит безобидно (1, -1/2, 1/6, -1/30…), то дальше эти числа начинают быстро — почти как факториалы! — расти (и для m=57 появляются коэффициенты ~10^{30})

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

***

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

про это поведение ряда Лейбница — взял из книжки «Mathematics by Experiment» (Borwein, Bailey) — и книжку вообще, и конкретно это место в ней мне показал Сережа Маркелов

а сама математика восходит к Эйлеру — которому очень хотелось посчитать 20 знаков суммы обратных квадратов, не складывая 10^{20} чисел

(после этого Эйлер нашел и точное значение этой суммы… а также ζ(2n), ζ(-k), и даже функциональное уравнение для ζ(x) — видимо поэтому ζ(x) теперь называют дзета-функцией… Римана)

сохраню небольшое видео с введением в Lean:

https://youtu.be/I2zaPoj3G50

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

Рассказывал пару раз на летних школах про количества решений разных уравнений mod p. Вот, например, на такую таблицу можно попробовать посмотреть и поискать какие-то закономерности.

Генерировал ответы в таком духе:

from prettytable import PrettyTable

def isPrime(n):
    if n == 2 or n == 3: return True
    if n%2 == 0 or n<2: return False
    for i in range(3, int(n**0.5)+1, 2):
        if n%i == 0: return False    
    return True

t = PrettyTable(['p','#{x²=-1}','#{x²=2}','#{x²+y²=1}',\
  '#{y²=x³+x²}','#{y²=x³-x}'])
for p in range(3,50+1,2):
 counts = [0]*5
 if isPrime(p):
            for x in range(p):
                if (x*x+1)%p == 0: counts[0] += 1
                if (x*x-2)%p == 0: counts[1] += 1
                for y in range(p):
                    if (x*x+y*y-1)%p == 0: counts[2] += 1
                    if (x*x*x+x*x-y*y)%p == 0: counts[3] += 1
                    if (x*x*x-x-y*y)%p == 0: counts[4] += 1
            t.add_row([p,*counts])
 else:
  t.add_row([p,'','','','',''])
t.align="r"
print(t)

Пример того, на что еще интересно смотреть (и при этом есть достаточно доступный ответ) — количества точек на сферах mod p, начиная с x²+y²+z²=1 mod p.

Если хочется менее расплывчатых вопросов — то вот такие, например, задачи давал: https://dev.mccme.ru/~merzon/leto2022/p-count-ex.pdf

Если, наоборот, хочется возвышенных разговоров — то вот (для школьников, не для экспертов) есть рассказа А.Кузнецова про теорему Ферма, гипотезу Римана и алгебраическую геометрию: https://youtu.be/ElFNGh4WF6Y

в конце поста https://news.1rj.ru/str/compmathweekly/6 был вопрос без ответа, можно ли так получить асимптотически правильную (пусть с точностью до константы) оценку на факториал

а сегодня коллега Устинов всё объяснил

// дальше техническое математическое, простите

напомню контекст: e^n = …+n^n/n!+…, поэтому e^n>n^n/n! (мы оценили сумму положительных чисел самым большим слагаемым), т.е. n!>(n/e)^n

так вот, если n^n/n!=S(n), то следующий член равен S(n)/(1+1/n), т.е. больше, чем S(n)e^{-1/n}; следующий равен S(n)/(1+1/n)(1+2/n), т.е. больше, чем S(n)e^{-1/n}e^{-2/n} и т.д. — если отойти на m шагов, то получится хотя бы S(n)e^{-m(m+1)/2n}

в частности, если отойти на (как и было предсказано :) на √n шагов, то получится хотя бы с⋅S(n) (для некоторого с, не зависящего от n) — и это дает оценку e^n > √n с S(n), т.е. n! > с√n (n/e)^n

можно доказать и оценку с другой стороны:
(1+1/n)(1+2/n)…(1+m/n)>1+1/n+…+m/n>1+m²/2n, поэтому если отойти на √n шагов, то слагаемое уменьшится хотя бы в полтора раза и т.д. — получается, что n! < C√n (n/e)^n

можно вытащить и полноценную формулу Стирлинга с использованием формулы Валлиса (оценивая центральный биномиальный коэффициент)

на Московской математической олимпиаде была сегодня такая задача

«У хозяйки есть кусок мяса, которым она хочет накормить трёх котиков. Раз в несколько секунд хозяйка отрезает кусочек мяса и скармливает его одному из котиков на свой выбор, причём каждый кусочек должен составлять одну и ту же долю куска, от которого его отрезают. Через некоторое время хозяйка убирает остаток мяса в холодильник. Может ли она скормить котикам поровну мяса?» (А.Кушнир)

другими словами, можно ли разбить члены какой-нибудь конечной геометрической прогрессии с 0<q<1 на 3 группы с одинаковыми суммами?

любая хозяйка, которая умеет резать мясо в золотом сечении, решит задачу для 2 котиков; для 3 котиков попробуйте решить на бумажке; для 4 котиков покажу попозже (пусть пока будет возможность подумать), как воспользоваться компьютером; какой ответ для 5 котиков — хотел бы знать

пусть мы хотим разбить геометрическую прогрессию 1,q,…,q^{N-1} на 4 равные части E, F, G, H

если E(q)=F(q), то многочлен F-E делится на минимальный многочлен числа q — и то же верно для G-E, H-E

но тогда делится на этот минимальнй многочлен и (F-E)+(G-E)+(H-E)=(E+F+G+H)-4E=1+q+…+q^{N-1}-4E

это условие совершенно не достаточное, только необходимое, но все же попробуем так протестировать гипотетическую долю для котика E=1:

import sympy
q = sympy.Symbol('q')
E = 1
for N in range(3,40):
    print(N,sympy.factor(-4*E+(1-q**N)/(1-q)))

видим, что нужный многочлен получается приводимым в обозримом диапазоне только дважды: для N=4 (этот случай явно не подходит) и для N=14, где он делится на q³+q²-1

теперь можно попробовать взять N=14, доли двух котиков 1 и q²+q³, а доли оставшихся двух котиков поискать прямолинейным перебором (раздаем оставшиеся куски всевозможными способами между G и H и проверяем, делятся ли G-1 и H-1 на q³+q²-1)

про это программу положу в комментарии

***

можно попробовать так же подойти к 5 котикам, но там для E=1 в обозримом диапазоне кандидатов на решение не видно!

можно пытаться перебирать разные E, но тут, наоборот, начинает лезть много паразитных решений (напр., 1+q^7, q+q^4, q^2+q^3, q^5+q^6 все равны при q=-1, но мясо это резать не помогает )

в общем, если сможете найти пример для деления на 5 частей, то пишите

***

впечатления от sympy… неоднозначные

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

нравится сюжет Конвея про аналогию между играми и числами

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

если в H выигрывает второй игрок, то результат у G+H такой же как и в G — это мотивирует объявить все выигрышные для второго игрока игры нулевыми

а вот игры, в которых выигрывает первый, бывают очень разными

если «ним-число» *n — это глуповатая игра «есть кучка из n камней, за ход можно взять любое количество камней из кучки», то *0 действительно нулевая игра, а все остальные *n — различные… и ненулевые )

и игра в четыре кучки камней *1+*3+*5+*7 уже не очень простая (не все персонажи фильма L'Année dernière à Marienbad справились), чтобы научиться в нее играть, хорошо бы изучить таблицу операций с ним-числами

вот такой, например, листок про это: https://dev.mccme.ru/~merzon/v14/pscache/5d-nim.pdf

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

def mex(N,arr):
    for a in range(N):
        if (a not in arr):
            return a
    return None

N = 2**(2**2)

t_sum = [list(range(N))]
for m in range(1,N):
    newline = []
    for i in range(N):
        # *m+*i = mex{*j+*i,*m+*i'|j<m,i'<i}
        arr = [line[i] for line in t_sum] + newline
        newline.append(mex(N,arr))
    t_sum.append(newline)
print(*t_sum,sep="\n")

t_mul = [[0]*N]
for m in range(1,N):
    newline = []
    for i in range(N):
        # *m.*i = mex{*j.(*i+*i')+*m.*i'|j<m,i'<i}
        arr = []
        for i1,mi1 in enumerate(newline):
            ii1 = t_sum[i][i1]
            for line in t_mul:
                jii1 = line[ii1] #*j.(*i+*i')
                arr.append(t_sum[jii1][mi1])
        newline.append(mex(N,arr))
    t_mul.append(newline)
print()
print(*t_mul,sep="\n")

можно заметить, а потом и доказать, что ним-сложение — это, на самом деле, просто побитовое сложение

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

( определение — можно прочитать в https://en.wikipedia.org/wiki/Nimber#Multiplication )

но операция оч. хорошая — в частности, ним-числа, меньшие *(2^(2^k)), образуют конечное поле

(мини-комментарий к предыдущему)

как расширить поле из 2 элементов до поля из 2^2=4 элементов? достаточно добавить элемент A, удовлетворяющий соотношению
A²=A+1

дальше получить поле из 2^(2^2) элементов можно расширением B²=B+A; потом можно взять расширение C²=C+B и т.д.

именно такую последовательность расширений и реализует ним-умножение (A=*(2^(2^0)), B=*(2^(2^1)), C=*(2^(2^2))…)

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