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
{
ld 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;
for (int a: A)
{
int d = 1;
for (int ai: A) if (a != ai) d *= a - ai;
res += (ld)(fast(a, b) - P.get(a)) / d;
}
P.x2 = res;
res = 0;
A.pop_back();
for (int a: A)
{
int d = 1;
for (int ai: A) if (a != ai) d *= a - ai;
res += (ld)(fast(a, b) - P.get(a)) / d;
}
P.x1 = res;
res = 0;
A.pop_back();
for (int a: A)
{
int d = 1;
for (int ai: A) if (a != ai) d *= a - ai;
res += (ld)(fast(a, b) - P.get(a)) / d;
}
P.x0 = res;
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.
