Всем привет. Недавно мне в голову пришла следующая задача: Рассмотрим множество точек (x, y) с целыми координатами таких, что 0 <= x < a и 0 <= y < b. Требуется найти количество остроугольных треугольников с вершинами в этих точках.
Попытки проинтерполировать
Ясно, что можно написать функцию f(a, b), которая будет искать ответ и работать при этом за (ab) ^ 3
. Я предположил, что она ведет себя как многочлен от двух переменных степени не более 6. Я попытался её проинтерполировать, используя эту теорему. Но у меня ничего не получилось, так как при мономах степени больше 6 интерполяция давала ненулевой коэффициент. Не получилось также с тупоугольными и прямоугольными треугольниками.
code (для прямоугольных треугольников)int stupid(int a, int b)
{
int ans = 0;
for (int x1=0;x1<a;x1++)
for (int x2=0;x2<a;x2++)
for (int x3=0;x3<a;x3++)
for (int y1=0;y1<b;y1++)
for (int y2=0;y2<b;y2++)
for (int y3=0;y3<b;y3++)
{
int a = (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2);
int b = (x1 - x3) * (x1 - x3) + (y1 - y3) * (y1 - y3);
int c = (x3 - x2) * (x3 - x2) + (y3 - y2) * (y3 - y2);
if ((a + b == c) && min(a, min(b, c))) ans++;
}
return ans / 2;
}
signed main()
{
int C1 = 5, C2 = 5;
ld res = 0;
for (int a1=4;a1<=C1+3;a1++)
{
for (int a2=4;a2<=C2+3;a2++)
{
int d = 1;
for (int ai=4;ai<=C1+3;ai++) if (ai != a1) d *= a1 - ai;
for (int ai=4;ai<=C2+3;ai++) if (ai != a2) d *= a2 - ai;
res += (ld)(stupid(a1, a2)) / d;
}
}
cout << setprecision(20) << fixed << res;
}
Данный код узнаёт коэффициент при мономе a ^ (C1 - 1) * b ^ (C2 - 1)
Что я хотел бы узнать:
- решается ли эта задача быстрее, чем за куб
- решается ли эта задача за O(1)
- может, кто-нибудь знает задачи, где формула для ответа не очевидна и её можно подобрать этим методом?
UPD: найдена формула для количества прямоугольных треугольников с b = 2: f(a, 2) = 2 * a ^ 2 - 4
, a > 1.
UPD2: большое спасибо bronze_coder за нахождение решения за O(1)
для b = const
: OEIS A189814.
Для интерполяции надо использовать ai > b ^ 2
. EDIT: ai > (b - 1) ^ 2
UPD3: Наконец, я написал решение за O((ab) ^ 2)
.
Codeint fast(int a, int b)
{
int ans = 0;
for (int a1=-a+1;a1<=a-1;a1++)
{
for (int b1=-b+1;b1<=b-1;b1++)
{
if (a1 == 0 && b1 == 0) continue;
for (int a2=-a+1;a2<=a-1;a2++)
{
for (int b2=-b+1;b2<=b-1;b2++)
{
if (a2 == 0 && b2 == 0) continue;
// первое условие - скалярное произведение равно 0
// второе условие - векторное произведение < 0 - чтобы не считать одну пару дважды
if (b1 * b2 + a1 * a2 == 0 && a1 * b2 - b1 * a2 < 0)
{
int cnta = a - max(max(abs(a1), abs(a2)), abs(a1 - a2));
int cntb = b - max(max(abs(b1), abs(b2)), abs(b1 - b2));
if (cnta > 0 && cntb > 0) ans += cnta * cntb;
}
}
}
}
}
return ans;
}
Теперь я могу использовать большие значения a
и b
для интерполяции.
Но всё равно кажется странным, что количество прямоугольных треугольников пропорционально (ab) ^ 2
, а не (ab) ^ 3
. Сейчас попробую понять, почему формула не работает для ai <= b ^ 2
. EDIT: ai <= (b - 1) ^ 2
UPD4: Код, который находит формулу для f(a, b) при a > b ^ 2
и работает за O(b ^ 6)
:
Спойлерint fast(int a, int b)
{
int ans = 0;
for (int a1=-a+1;a1<=a-1;a1++)
for (int b1=-b+1;b1<=b-1;b1++)
{
if (a1 == 0 && b1 == 0) continue;
for (int a2=-a+1;a2<=a-1;a2++)
for (int b2=-b+1;b2<=b-1;b2++)
{
if (a2 == 0 && b2 == 0) continue;
if (b1 * b2 + a1 * a2 == 0 && a1 * b2 - b1 * a2 < 0)
{
int cnta = a - max(max(abs(a1), abs(a2)), abs(a1 - a2));
int cntb = b - max(max(abs(b1), abs(b2)), abs(b1 - b2));
if (cnta > 0 && cntb > 0) ans += cnta * cntb;
}
}
}
return ans;
}
struct polynom
{
int x0 = 0, x1 = 0, x2 = 0;
void print()
{
cout << x2 << "x^2 + " << x1 << "x - " << -x0 << endl;
}
ld get(ld x)
{
return x2 * x * x + x1 * x + x0;
}
};
void interpolate(int b)
{
polynom P = {0, 0, 0};
vector<int> A = {b * b + 1, b * b + 2, b * b + 3};
ld res = 0;
vector<ld> B = {fast(A[0], b), fast(A[1], b), fast(A[2], b)};
for (int i=0;i<3;i++)
{
int d = 1;
for (int j=0;j<3;j++) if (A[i] != A[j]) d *= A[i] - A[j];
res += (B[i] - P.get(A[i])) / d;
}
P.x2 = res;
res = 0;
A.pop_back();
for (int i=0;i<2;i++)
{
int d = 1;
for (int j=0;j<2;j++) if (A[i] != A[j]) d *= A[i] - A[j];
res += (B[i] - P.get(A[i])) / d;
}
P.x1 = res;
P.x0 = B[0] - P.get(A[0]);
P.print();
}
Пытаясь найти закономерность в коэффициентах многочлена P, я обратился к OEIS, но ничего там не нашел :(, кроме x2 = b * (b - 1)
, что и так было очевидно.
UPD5:
Наконец-то нашёл формулу и решение за O(min(a, b) ^ 6)
Если a < b
, то поменяем их местами. Теперь нужно решить за O(b ^ 6)
. Если a <= (b - 1) ^ 2 + 1
, то запустим решение за O((ab) ^ 3)
= O(b ^ 6)
. Теперь нужно разобраться с большими a.
Определения
Рассмотрим все прямоугольные треугольники и разделим их на 3 типа:
треугольники, у которых нет вершин на прямой x = a - 1
треугольники, у которых есть вершины на прямой x = a - 1
и две стороны параллельны осям координат
остальные треугольники
Обозначим количество треугольников третьего типа за C(a)
(какая-то функция от a).
Количество треугольников второго типа можно посчитать по формуле (a - 1) * b * (b - 1) / 2 * 4
:
ДоказательствоВыберем одну вертикальную прямую x = z, z < a - 1
(a - 1
способ), вертикальную прямую x = a - 1
и две горизонтальные прямые (b * (b - 1) / 2
способов). Пересечением наших прямых являются 4 точки. Мы можем выкинуть из них любую, и после этого оставшиеся 3 точки будут образовывать прямоугольный треугольник второго типа. Ясно, что любой треугольник второго типа может быть представлен в виде пересечения таких прямых и отбрасывания одной из точек пересечения.
Из определения следует, что для вершин (x, y) треугольников первого типа выполняется 0 <= x < a - 1
и 0 <= y < b
, то есть их количество равно f(a - 1, b)
, по определению функции f
.
Итак, f(a, b) = f(a - 1, b) + (a - 1) * b * (b - 1) / 2 * 4 + C(a)
Теперь докажем, что C(a)
при всех a > (b - 1) ^ 2 + 1
— константа.
ДоказательствоЛемма. Максимальный модуль разности абсцисс вершин треугольника третьего типа не превосходит (b - 1) ^ 2 + 1
Для доказательства посмотрим на рисунок.
Прикрепился только повёрнутый( У нашего треугольника есть хотя бы две стороны, непараллельные осям координат. Пусть они AB и BC. Пусть угловой коэффициент прямой AB равен k, тогда у прямой BC он равен -1/k
. Тогда разность абсцисс верши A и C не больше, чем (b - 1) / |k| + (b - 1) / |-1/k| = (b - 1) * (|k| + 1 / |k|)
.
Известно, что если мы возьмём два положительных числа и будем их сближать на числовой прямой с сохранением произведения, то их сумма будет уменьшаться. Тогда для максимизации |k| + 1 / |k|
нужно минимизировать или максимизировать |k|. Из ограничения по y-координатам, |k| < b - 1
и 1 / |k| < b - 1
. Тогда максимальное значение нашего выражения (b - 1) * ((b - 1) + 1 / (b - 1)) = (b - 1) ^ 2 + 1
. Лемма доказана.
Нужно доказать, что при a > (b - 1) ^ 2 + 1
количество треугольников третьего типа не зависит от a. Это очевидно следует из леммы, ведь каждому такому треугольнику для a
мы можем сопоставить треугольник для a + 1
, полученный из первого сдвигом на 1 вправо. И каждому треугольнику 3 типа для a + 1 аналогично сопоставим треугольник 3 типа для a, сдвинув его на 1 влево.
C(a)
— константа. Обозначим эту константу за c
.
Итак, f(a, b) = f(a - 1, b) + (a - 1) * b * (b - 1) / 2 * 4 + c
. Немного преобразований, и мы получаем формулу для f(a, b)
через f((b - 1) ^ 2 + 1, b)
.
Реализацияll fast(int a, int b)
{
ll ans = 0;
for (int a1=-a+1;a1<=a-1;a1++)
for (int b1=-b+1;b1<=b-1;b1++)
{
if (a1 == 0 && b1 == 0) continue;
for (int a2=-a+1;a2<=a-1;a2++)
for (int b2=-b+1;b2<=b-1;b2++)
{
if (a2 == 0 && b2 == 0) continue;
if (b1 * b2 + a1 * a2 == 0 && a1 * b2 - b1 * a2 < 0)
{
int cnta = a - max(max(abs(a1), abs(a2)), abs(a1 - a2));
int cntb = b - max(max(abs(b1), abs(b2)), abs(b1 - b2));
if (cnta > 0 && cntb > 0) ans += cnta * cntb;
}
}
}
return ans;
}
struct polynom
{
ll x0 = 0, x1 = 0, x2 = 0;
void print()
{
cout << x2 << "x^2 + " << x1 << "x - " << -x0 << endl;
}
ld get(ld x)
{
return x2 * x * x + x1 * x + x0;
}
};
ll get(int a, int b)
{
if (a < b) swap(a, b);
if (a <= (b - 1) * (b - 1) + 1) return fast(a, b);
int x = fast((b - 1) * (b - 1) + 1, b);
int y = fast((b - 1) * (b - 1) + 2, b);
// y = x + c + (a - 1) * b * (b - 1) / 2 * 4;
int c = y - x - ((b - 1) * (b - 1) + 1) * (b - 1) * b * 2;
polynom P;
P.x2 = b * b - b;
P.x1 = c - b * b + b;
P.x0 =-b*b*b*b*b*b + 5*b*b*b*b*b - 11*b*b*b*b + 13*b*b*b + b*b*(-8 - c) + b * (2 * c + 2) - 2 * c + x;
return P.get(a);
}
К сожалению, c
для разных b
— разные, и я не смог найти между ними закономерность. Не помогли ни интерполяция, ни OEIS. Осталось несколько вещей, которые надо сделать:
- Выразить
c
через b
, ведь c
зависит только от b
- Решить задачу для
a <= (b - 1) ^ 2
быстрее - Решить задачу для остроугольных треугольников
Будет жёстко.
UPD6: Можно ускорить решение до O(b^5)
, используя эту идею
Auto comment: topic has been translated by polosatic (original revision, translated revision, compare)
The right triangle problem is OEIS A077435
EDIT: OEIS A189814 gives the answer for rectangles.
Thank you! I was quite sure that problem with right triangles can be solved in
O(a * b)
But unfortunately, f(n, n) also is not a polynomial of one variable.OMG! I used too small
ai
for interpolation with constant b, so it didn't work. After usingai >= b ^ 2
everything works!Можно как минимум за квадрат решить, перебирая все пары направляющих векторов.
Действительно! Как я не догадался.
Автокомментарий: текст был обновлен пользователем polosatic (предыдущая версия, новая версия, сравнить).
Very nice
A not-implemented solution: let the triangle be tree vectors $$$a,b,c$$$.brute force $$$b-a$$$, and you will get $$$\dfrac{c-a}{|c-a|}$$$(which is rotate it for 90 deg(you can asumme ccw or cw))(In program, you can do it by a ). then you can brute force $$$c-a$$$ by $$$O(a + b)$$$. Then the problem is:
which is trival.
If I understand you correctly, I have already implemented this solution in UPD3 of this blog
No. This solution is $$$O(ab(a+b))$$$.
UPD3 is
Brute force b-a and c-a, find if b-a ⊥ c-a and find how many triangles with vector b-a and c-a.
But we can just brute force $$$|c-a|$$$(The length of c-a). the angle of $$$c-a$$$ is found by rotate $$$b-a$$$.
Wow! Good idea! Thank you.
why $$$O(a^2b^2(a+b))$$$ is $$$O(b^5)$$$? may $$$a \gg b^2$$$.
The solution in UPD5 needs to calculate
f((b - 1) ^ 2 + 1, b)
andf((b - 1) ^ 2 + 2, b)
to calculatef(a, b)
, it works inO(b ^ 6)
, but can be optimized toO(b ^ 5)
using your method.