Точный подсчёт прямоугольных треугольников с вершинами в узлах сетки n×n: алгоритмы со сложностью O(n^(3/2+ε)) и O(n log² n)

Точный подсчёт прямоугольных треугольников с вершинами в узлах сетки n×n: алгоритмы со сложностью O(n^3/2+ε) и O(n log² n).

Аннотация. Рассматривается задача точного вычисления числа R(n) прямоугольных треугольников с вершинами в узлах сетки n×n (последовательность A077435 в OEIS). Получена формула в виде четырёхкратной суммы R(n) = 4S² + 4Σ (n − gp − μq)₊(n − max(gq, μp))₊ по p, q, g, μ ≥ 1 при gcd(p, q) = 1. Слагаемое симметрично при обмене пар (g, μ) и (p, q); из этого выводится алгоритм со сложностью O(n^3/2+ε). Далее условие взаимной простоты снимается функцией Мёбиуса, сумма разбивается на слои, каждый слой покрывается двумя малыми парами переменных, а оставшиеся двумерные суммы по многоугольникам вычисляются ядром сумм целых частей с моментами за O(log n); общая сложность составляет O(n log² n) арифметических операций. Обе реализации на Python совпали с прямым перебором для n ≤ 7, со всеми значениями OEIS для n ≤ 35, с b-файлом OEIS на 161 значении до n = 10 000 (алгоритм O(n^3/2+ε)) и на 69 значениях до n = 3 000 (почти линейный алгоритм). В измерениях почти линейный алгоритм оказался быстрее при n = 12 800 (9,35 с против 11,55 с).

1. Введение

Число R(n) прямоугольных треугольников, вершины которых лежат в узлах решётки {1, ..., n}×{1, ..., n}, составляет последовательность A077435 [1]. В OEIS приведены значения до n = 10 000, вычисленные с помощью таблицы чисел треугольников с заданным ограничивающим прямоугольником (последовательность A279433 [2]); закрытой формулы для R(n) в OEIS нет.

Цель работы - вывести точную формулу, построить два алгоритма вычисления R(n) для одного n и экспериментально измерить число операций и время. Все вычисления выполняются в целых числах.

2. Постановка задачи

Пусть n ≥ 1 и S = n(n − 1)/2. Требуется вычислить R(n) - число неупорядоченных троек узлов {1, ..., n}², являющихся вершинами прямоугольного треугольника. В работе сравниваются три реализации:

базовая, в которой вычисляется четырёхкратная сумма (1) прямым перебором всех допустимых четвёрок (p, q, g, μ); число слагаемых порядка n²log² n;
алгоритм O(n^3/2+ε), использующий симметрию слагаемого (раздел 4);
почти линейный алгоритм O(n log² n), использующий слои Мёбиуса и ядро сумм целых частей (раздел 5).

Все три реализации возвращают одно и то же целое число; эксперимент направлен на оценку влияния способа суммирования.

3. Формула

Теорема 1. Для всех n ≥ 1

R(n) = 4S² + 4W, W = Σ_{p,q≥1, gcd(p,q)=1} Σ_g,μ≥1 (n − gp − μq)₊·(n − max(gq, μp))₊,

(1)

где x₊ = max(x, 0).

Доказательство. Фиксируем вершину B прямого угла и катеты u = A − B, v = C − B (ненулевые, u ⊥ v). Каждый треугольник соответствует двум упорядоченным парам (u, v). Тройка B, B+u, B+v помещается в сетку числом способов (n − X)₊(n − Y)₊, где X×Y - размеры ограничивающего прямоугольника точек 0, u, v. Поэтому R(n) = (1/2)Σ (n − X)₊(n − Y)₊ по упорядоченным перпендикулярным парам.

Любой ненулевой вектор записывается как u = g·(p, q), где (p, q) - примитивный вектор, g ≥ 1; перпендикулярный ему вектор имеет вид v = λ·(−q, p), λ ≠ 0. Повороты на 90° и отражения сетки сохраняют вес (n − X)(n − Y), поэтому сумма равна 4-кратной сумме по примитивным векторам с p ≥ 1, q ≥ 0. Для p, q ≥ 1 при λ > 0 получаем X = gp + λq, Y = max(gq, λp); при λ = −μ < 0 получаем X = max(gp, μq), Y = gq + μp, и второй случай переходит в первый обменом p ↔ q. Направления (1, 0) и (0, 1) дают катеты вдоль осей и вклад 2S² в сумму, то есть слагаемое 4S² в R(n). Собирая вклады и сокращая множитель 1/2, получаем (1).

4. Симметрия и алгоритм O(n^3/2+ε)

Лемма 1. Слагаемое (n − gp − μq)(n − max(gq, μp)) не меняется при обмене пар (g, μ) ↔ (p, q).

Доказательство. При замене g ↔ p, μ ↔ q выражение gp + μq переходит в pg + qμ, а набор {gq, μp} - в {pμ, qg}, то есть в тот же набор.

Лемма 2. Пусть T = ⌊√n⌋. Для ненулевых слагаемых выполнено gp, gq, μp, μq < n; поэтому либо g, μ < T, либо p, q ≤ (n − 1)//T.

Доказательство. Если g ≥ T или μ ≥ T, то из gp, gq < n (соответственно μp, μq < n) получаем p, q ≤ (n − 1)//T.

Алгоритм. Случай I (g, μ < T): перебираются g, μ и p, а сумма по q, взаимно простым с p, вычисляется по формулам для степенных сумм с мёбиусовской суммой по свободным от квадратов делителям числа p; точки разделения q по условию gq > μp вычисляются непосредственно. Случай II (g ≥ T или μ ≥ T): перебираются взаимно простые p, q ≤ (n − 1)//T и g, а сумма по μ вычисляется по формулам. Случаи не пересекаются.

Теорема 2. Описанный алгоритм вычисляет R(n) за O(n^3/2+ε) арифметических операций.

Доказательство. Число узлов перебора в случае I не превосходит Σ_g,μ<T n/max(g, μ) ≤ 2n^3/2, в случае II - не более 2n^3/2(1 + o(1)). Стоимость узла в случае I не превосходит числа свободных от квадратов делителей p, то есть 2^ω(p) = n^o(1), в случае II она постоянна.

5. Почти линейный алгоритм

Слои Мёбиуса. Заменяя условие gcd(p, q) = 1 мёбиусовской суммой по d = gcd, получаем W = Σ_d≥1 μ(d) S_d, где S_d = Σ_g,μ,p,q≥1 (n − dX)₊(n − dY)₊, X = gp + μq, Y = max(gq, μp) (без условия взаимной простоты). Ненулевые слагаемые удовлетворяют X, Y ≤ N_d = ⌊(n − 1)/d⌋, а вес (n − dX)(n − dY) = n² − nd(X + Y) + d²XY. Следовательно,

S_d = n²Z₀(N_d) − nd·Z₁(N_d) + d²Z₂(N_d),

(2)

где Z₀(N), Z₁(N), Z₂(N) - суммы величин 1, X + Y, XY по всем четвёркам (g, μ, p, q) ≥ 1 с X, Y ≤ N. Значение N_d принимает O(√n) различных значений; d с одинаковым N_d объединяются с помощью префиксных сумм Σμ(d), Σdμ(d), Σd²μ(d), вычисляемых решетом за O(n).

Покрытие двумя малыми парами. Пусть B = ⌊√N⌋. Если g ≥ B + 1 или μ ≥ B + 1, то p, q ≤ B (как в лемме 2). Обозначим G_k(g, μ; N) сумму по (p, q) ≥ 1 из области X, Y ≤ N, а G_k^(B) - ту же сумму с дополнительными условиями p, q ≤ B. По лемме 1

Z_k(N) = 2Σ_g,μ=1^B G_k(g, μ; N) − Σ_g,μ=1^B G_k^(B)(g, μ; N).

(3)

Действительно, множество всех четвёрок есть объединение {g, μ ≤ B} и {g > B или μ > B}; вторая часть лежит в {p, q ≤ B}, и по лемме 1 её вклад равен вкладу {g, μ ≤ B} минус вклад четвёрок, у которых все четыре переменные не превосходят B. Дополнительно G(g, μ) = G(μ, g) (обмен g ↔ μ и p ↔ q не меняет X и Y), поэтому достаточно перебрать пары g ≤ μ. Число таких пар равно B(B + 1)/2.

Многоугольные суммы. Для фиксированных g, μ область суммирования по (p, q) - многоугольник, ограниченный прямыми p = 1, q = 1, gp + μq = N и ограничениями μp ≤ N, gq ≤ N; линия gq = μp разбивает его на две части с Y = μp и Y = gq. Суммирование по q даёт степенные суммы с границами вида ⌊(ap + b)/m⌋; остаётся сумма по p выражений pⁱ⌊(ap + b)/m⌋^j, i + j ≤ 3. Точки излома границ находятся бисекцией по монотонным предикатам за O(log N). Суммы pⁱ⌊(ap + b)/m⌋^j с шестью состояниями (i, j) вычисляются рекурсией Евклида за O(log N); формулы взаимного шага выведены из тождества Σ_x xⁱf(x)^j = Σ_t ((t+1)^j − t^j)·(P_i(n) − P_i(g(t)+1)).

Теорема 3. Алгоритм вычисляет R(n) за O(n log² n) арифметических операций.

Доказательство. На слой N затрачивается B(B + 1)/2 ≤ N пар, на каждую пару - O(1) вызовов ядра по O(log n) операций. Различных значений N вида ⌊(n − 1)/d⌋ существует O(√n): при d ≤ √n суммарное число пар не превосходит Σ_d≤√n n/(2d) = (n/4)ln n + O(n), а значения N ≤ √n в сумме дают O(n). Итого O(n log n) пар и O(n log² n) операций; решето для μ(d) требует ещё O(n log log n).

6. Методика эксперимента и проверка корректности

Все реализации написаны на Python и используют только стандартную библиотеку; целые числа произвольной длины. Время измерялось функцией time.perf_counter; для n ≤ 800 приведён минимум из трёх повторов, для больших n - один запуск. Число операций: для базовой реализации - число обработанных четвёрок, для алгоритма O(n^3/2+ε) - число узлов перебора, для почти линейного - число пар (g, μ) ≤ B, для которых вычислялись многоугольные суммы.

Корректность проверялась так:

для n ≤ 7 все реализации совпали с прямым перебором троек узлов и с OEIS;
для n ≤ 35 алгоритм O(n^3/2+ε) и почти линейный алгоритм совпали со всеми значениями OEIS, приведёнными на странице A077435;
для всех n ≤ 150 три реализации совпали между собой;
алгоритм O(n^3/2+ε) совпал с b-файлом OEIS на 161 значении (в том числе n = 100, 200, 300, 500, 800, 1 000, 2 000, 3 000, 5 000, 8 000, 10 000 и 150 случайных n ≤ 1 500), расхождений нет;
почти линейный алгоритм совпал с b-файлом OEIS на 69 значениях до n = 3 000, расхождений нет;
в измерениях раздела 7 значения всех реализаций, запускавшихся для данного n, проверялись на равенство.

Вспомогательные части (ядро сумм целых частей, многоугольные суммы G, покрытие (3)) отдельно проверялись перебором на тысячах случайных наборов параметров.

7. Результаты

n	Базовый метод, с	Четвёрок	Метод n^3/2, с	Узлов	Почти линейный метод, с	Пар
50	0,002	11 782	0,001	760	0,010	66
100	0,011	66 239	0,003	2 361	0,025	135
200	0,064	353 587	0,009	7 201	0,064	311
400	0,384	1 812 003	0,034	21 432	0,175	652
800	1,782	9 026 069	0,097	62 636	0,323	1 386
1 600	—	—	0,308	182 568	0,768	2 972
3 200	—	—	0,943	525 619	1,801	6 257
6 400	—	—	2,969	1 510 598	4,031	13 161
12 800	—	—	11,554	4 318 097	9,346	27 798

Для n = 800 алгоритм O(n^3/2+ε) быстрее базового в 18,4 раза (0,097 с против 1,782 с) при сокращении числа слагаемых с 9 026 069 до 62 636. При n ≤ 6 400 алгоритм O(n^3/2+ε) быстрее почти линейного; при n = 12 800 почти линейный алгоритм впервые оказывается быстрее: 9,346 с против 11,554 с (в 1,24 раза). Оценка точки пересечения по аппроксимации степенными законами даёт n порядка 9 800.

Эмпирические показатели степени, найденные методом наименьших квадратов в логарифмическом масштабе: время базового метода - 2,45 (n = 100...800), число его слагаемых - 2,36; время алгоритма O(n^3/2+ε) - 1,73 (n = 1 600...12 800), число его узлов - 1,52; время почти линейного алгоритма - 1,20, число пар - 1,08. Число пар для больших n вычислялось аналитически:

n	Число пар	Пар / (n ln n)
10³	2 255	0,3264
10⁵	341 702	0,2968
10⁷	45 685 958	0,2834
10⁹	5 719 404 705	0,2760

Отношение монотонно убывает к значению 1/4, предсказанному оценкой (n/4)ln n из теоремы 3; показатель, найденный по аналитическим данным (n = 10⁵...10⁹), равен 1,06. Время алгоритма O(n^3/2+ε) растёт быстрее числа его узлов (показатель 1,73 против 1,52), так как при больших n возрастает стоимость арифметики с большими целыми и работы с кэшем делителей. Значение R(12800) = 201 480 946 727 162 368 получено обоими алгоритмами и лежит за пределами b-файла OEIS (n ≤ 10 000).

8. Заключение

В работе получена формула (1) для числа прямоугольных треугольников в сетке n×n в виде четырёхкратной суммы и два алгоритма вычисления R(n) для одного n. Алгоритм O(n^3/2+ε) основан на симметрии слагаемого относительно обмена пар (g, μ) ↔ (p, q). Почти линейный алгоритм O(n log² n) использует слои Мёбиуса, покрытие двумя малыми парами и ядро сумм целых частей.

Обе реализации совпали с b-файлом OEIS в проверенных точках (161 значение до n = 10 000 и 69 значений до n = 3 000 соответственно). На Python для n = 800 алгоритм O(n^3/2+ε) быстрее базового в 18,4 раза, а при n = 12 800 почти линейный алгоритм быстрее алгоритма O(n^3/2+ε) в 1,24 раза.

Исходные тексты реализаций, тесты корректности и скрипт измерений доступны в репозитории: https://github.com/tankwars92/right_triangles.

Список использованных источников

1. The On-Line Encyclopedia of Integer Sequences. A077435: Number of right triangles whose vertices are lattice points in an n X n grid of points [Электронный ресурс]. - URL: https://oeis.org/A077435 (дата обращения: 17.10.2026).

2. The On-Line Encyclopedia of Integer Sequences. A279433 [Электронный ресурс]. - URL: https://oeis.org/A279433 (дата обращения: 17.10.2026).