Hello, Codeforces. Recently I invented a problem. Consider a set of points (x, y) with integer coordinates, 0 <= x < a и 0 <= y < b. We want to know the number of acute triangles with vertices at given points.
Trying to interpolate
We can write function f(a, b) which finds the answer and works in (ab) ^ 3
. I suggested that it is a polynomial of two variables with degree <= 6. I tried to interpolate it using this theorem. But my code finds non-zero coefficients with monoms with degrees > 6, so my hypothesis was not confirmed. I also failed with right triangles.
code (for right triangles)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;
}
This code finds a coefficient with a ^ (C1 - 1) * b ^ (C2 - 1)
What I would like to know:
- can this problem be solved faster than O(stupid(a, b))
- can this problem be solved in O(1)
- maybe, someone knows problems with difficult to find formulas and this method can help to find them?
UPD: Found formula for the number of right triangles with b = 2: f(a, 2) = 2 * a ^ 2 - 4
, a > 1.
UPD2: Many thanks to bronze_coder for finding O(1)
solution for constant b: OEIS A189814.
Interpolation should use ai > b ^ 2
.
UPD3: Finally wrote a O((ab) ^ 2)
solution.
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;
// first condition is that dot product = 0
// second condition is that cross product < 0 not to count the same pair twice
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;
}
Now I can interpolate using bigger values of a
and b
.
But it is quite strange for me that the number of right triangles is directly proportional to (ab) ^ 2
instead of (ab) ^ 3
. Now I will try to understand why the formula doesn't work with ai <= b ^ 2
.
UPD4: Code that finds a formula for f(a, b) for a > b ^ 2
and works in O(b ^ 6)
:
Spoilerint 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();
}
I used OEIS to find a regularity between the coefficients of P, but i didn't find anything :( except x2 = b * (b - 1)
, but it was obvious.
UPD5:
Finally found a formula and solution in O(min(a, b) ^ 6)
If a < b
, let's swap them. Now we need to solve it in O(b ^ 6)
. If a <= (b - 1) ^ 2 + 1
, we can run solution in O((ab) ^ 3)
= O(b ^ 6)
. Now we need to solve it for a > (b - 1) ^ 2 + 1
.
Definitions
Consider all right triangles and divide them into 3 types:
triangles without any vertices on the line x = a - 1
triangles with some vertices on the line x = a - 1
and with two sides which are parallel to the x-axis or y-axis
other triangles
Let C(a)
be the number of third type triangles (some function).
The number of second type triangles is (a - 1) * b * (b - 1) / 2 * 4
:
ProofLet's choose a vertical line x = z, z < a - 1
(a - 1
variants), a vertical line x = a - 1
and two horizontal lines (b * (b - 1) / 2
variants). Intersection of these lines consists of 4 points. We can throw away any of them and the triangle formed from other 3 points has the second type. Obviously, every second type triangle can be presented as an intersection of these lines without one of points.
By definition, for every vertex (x, y) of first type triangles, 0 <= x < a - 1
and 0 <= y < b
, then the number of the triangles is f(a - 1, b)
, by definition of f
.
Well, f(a, b) = f(a - 1, b) + (a - 1) * b * (b - 1) / 2 * 4 + C(a)
Let's prove that C(a)
for every a > (b - 1) ^ 2 + 1
has same values.
ProofLemma. Maximal absolute value of difference between abscissas of two vertices of a third type triangle does not exceed (b - 1) ^ 2 + 1
Let's look on a picture
I could add it only rotated( The triangle has at least two sides which are parallel to neither x-axis nor y-axis (AB and BC). Let k be the slope of line AB, then the slope of BC is equal to -1/k
. Then, the difference between abscissas of A and C does not exceed
(b - 1) / |k| + (b - 1) / |-1/k| = (b - 1) * (|k| + 1 / |k|)
.
It is a well-known fact that if we have two positive numbers and we will bring their values together and their product remains the same, then their sum will decrease. Then, to maximize |k| + 1 / |k|
, we need minimize or maximize |k|. According to limitation of coordinates, |k| < b - 1
and 1 / |k| < b - 1
. Then the maximum value is (b - 1) * ((b - 1) + 1 / (b - 1)) = (b - 1) ^ 2 + 1
. The Lemma is proved.
We need to prove that for a > (b - 1) ^ 2 + 1
the number of third type triangles does not depend on a. It is obvious after proving Lemma because we can compare every triangle for a
with a triangle for a + 1
by shifting it to right.
C(a)
is a constant. Let c
be this constant.
Well, f(a, b) = f(a - 1, b) + (a - 1) * b * (b - 1) / 2 * 4 + c
. After some transforms, we get a formula for f(a, b)
using f((b - 1) ^ 2 + 1, b)
.
Implementationll 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);
}
Unfortunately, the value of c
is different for different values of b
and I could not find a regularity between them. Interpolation and OEIS did not help me. There are some things I should do:
- Express
c
through b
, after all c
depends on b
only - Solve the problem for
a <= (b - 1) ^ 2
faster - Solve the problem for acute triangles.
It will be hard.
UPD6: We can speed up the solution to O(b^5)
using this idea
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!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.