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

Как известно, сфера с g ручками получается из 4g-угольника попарной склейкой сторон. Возьмем (при g >= 2) 4g-угольник в плоскости Лобачевского, у которого отождествленные пары сторон равны по длине, а сумма углов равна 2π. Тогда на склеенной из него поверхности остается гиперболическая метрика; обратно, любую гиперболическую поверхность можно разрезать (по геодезическим отрезкам) на 4g-угольник, вершины которого сойдутся в выбранной на поверхности точке.

Если теперь мы рассмотрим универсальное накрытие поверхности плоскостью Лобачевского, то соответствующий 4g-угольник замостит плоскость своими копиями. Так мы можем изображать гиперболические метрики на поверхности (точки в пространстве Тейхмюллера) при помощи замощений.

Я написал небольшую программу, которая рисует фрагмент случайного такого замощения для g=2, т.е. копиями восьмиугольника. Склеиваемые стороны обозначены одинаковым цветом.

Я существенно опирался на архивные материалы аналогичного (но имевшего больший интерактивный функционал) проекта "The Teichmüller Navigator", созданного в 1993 Deva Van Der Werf под руководством David Ben-Zvi и Paul Burchard - каким-то чудом страница проекта сохранилась до сих пор: http://www.geom.uiuc.edu/apps/teich-nav/report/report.html.

ссылка на github с исходным кодом и инструкциями: https://github.com/mntroshkin/teichmuller-explorer.

надеюсь, что «гостевые посты» коллег здесь — тоже станут традицией

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

обычно ничего особо не видно, но прочитал недавно у Проппа, что поучительно рассматривать разбиения квадрата 2N×2N на доминошки и квадратики 2×2 по модулю степеней двойки

и действительно, компутер быстро посчитал мне начало последовательности
1
3
165
205879
5811552169
3676802302990923
51863128477992833745597 16264377685372072288158612804415

выглядит устрашающе, но по модулю 8, скажем, это просто
1, 3, 5, 7, 1, 3, 5, 7…

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

(технически модифицировал чуть-чуть как раз код для подсчета разбиений на доминошки динамикой по профилю из https://www.group-telegram.com/compmathweekly.com/17 — и немножко лениво сделал, мб из-за этого не очень быстро работает)

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

один из источников — статья Проппа «Some 2-adic conjectures concerning polyomino tilings of Aztec diamonds» (https://arxiv.org/abs/2204.00158)

работа над ошибками

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

import flint

def count_tilings(n,m,ext):
    ans = flint.fmpz_mat(2**n,1)
    ans[0,0] = 1
    for _ in range(m):
        ans = ext * ans
    return ans[0,0]

def newext(n,ext):
    ans = flint.fmpz_mat(2**n,2**n)
    for mask in range(2**n):
        for perp in range(2**n):
            if mask&perp == 0:
                if (mask+perp)%2 == 1:
                    ans[mask,perp] = ext[n-1][mask>>1,perp>>1]
                elif (mask+perp)%4 == 0:
                    ans[mask,perp] = ext[n-2][mask>>2,perp>>2]
    return ans

N = 5
ext = [flint.fmpz_mat(1,1,[1]),flint.fmpz_mat(2,2,[0,1,1,0])]
for n in range(2,2*N+1):
    ext.append(newext(n,ext))
    if n%2 == 0:
        ans = count_tilings(n,n,ext[n])
        print(n,ans)

12988816 разбиений квадрата 8×8 теперь подсчитываются за пару сотых секунды



теперь легко проверить, что разные делимости есть и для разбиений прямоугольников на доминошки — например, для 6×2 разбиений 13 => для (7k+6)×(3l+2) количество разбиений делится на 13 (и действительно, для 13×20 получается 5422065916132126528319352874496 разбиений, это число делится на 13)

количество разбиений прямоугольника 2×N — это просто число Фибоначчи, и что F_N | F_{2N+1} легко доказать комбинаторно (посмотреть как там покрыт доминошками средний столбец)

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

***

еще переписал генерацию случайного разбиения на доминошки ацтекского брильянта — через domino shuffling (aka urban renewal)

полный код не влезает в лимит телеграма, так что положу только основную функцию, которая делает из случайного брильянта порядка size случайный брильянт порядка size+1… point, наверное, в том, что получилось не сложнее, чем наивный код с флипами (но, конечно, флипы работают для разных областей, а не только для брильянтов)

def extend(field,size):
    newfield = [ [' ']*(2*size+2) for _ in range(2*size+2) ]
    for i in range(2*size+2):
        for j in range(2*size+2):
            if inside(i,j,size+1) and newfield[i][j]==" ":
                if inside(i,j-1,size) and field[i][j-1] in "Uu" and (not inside(i-1,j-1,size) or not field[i-1][j-1] in "Dd"):
                    newfield[i][j] = field[i][j-1]
                elif inside(i-2,j-1,size) and field[i-2][j-1] in "Dd" and (not inside(i-1,j-1,size) or not field[i-1][j-1] in "Uu"):
                    newfield[i][j] = field[i-2][j-1]
                elif inside(i-1,j,size) and field[i-1][j] in "Ll" and (not inside(i-1,j-1,size) or not field[i-1][j-1] in "Rr"):
                    newfield[i][j] = field[i-1][j] 
                elif inside(i-1,j-2,size) and field[i-1][j-2] in "Rr" and (not inside(i-1,j-1,size) or not field[i-1][j-1] in "Ll"):
                    newfield[i][j] = field[i-1][j-2]
                else:
                    bit = getrandbits(1)
                    newfield[i][j], newfield[i][j+1] = ("U", "u") if bit else ("L", "R")
                    newfield[i+1][j], newfield[i+1][j+1] = ("D", "d") if bit else ("l", "r")
    return newfield

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