Мои попытки ускорить куб до O(1)

#	User	Rating
1	tourist	3690
2	jiangly	3647
3	Benq	3581
4	orzdevinwang	3570
5	Geothermal	3569
5	cnnfls_csy	3569
7	Radewoosh	3509
8	ecnerwala	3486
9	jqdai0815	3474
10	gyh20	3447

#	User	Contrib.
1	maomao90	174
2	awoo	165
3	adamant	161
4	TheScrasse	160
5	nor	158
6	maroonrk	156
7	-is-this-fft-	152
8	orz	146
8	SecondThread	146
10	pajenegod	145

Всем привет. Недавно мне в голову пришла следующая задача: Рассмотрим множество точек (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).

Code

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;

		    // первое условие - скалярное произведение равно 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, 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), используя эту идею

Rev.	By	When	Δ	Comment
en16	polosatic	2023-02-07 09:55:02	82
ru18	polosatic	2023-02-07 09:54:01	93
ru17	polosatic	2023-02-07 05:28:56	129
en15	polosatic	2023-02-07 05:27:56	131
en14	polosatic	2023-02-07 05:25:34	35	Tiny change: 's together, then the' -> 's together and their product remains the same, then the'
ru16	polosatic	2023-02-06 19:10:33	6
en13	polosatic	2023-02-06 19:09:11	4916	Tiny change: ' < b`, let`s swap the' -> ' < b`, let's swap the'
ru15	polosatic	2023-02-06 18:23:23	26	Мелкая правка: 'реугольник. Ясно, чт' -> 'реугольник второго типа. Ясно, чт'
ru14	polosatic	2023-02-06 17:59:05	5245	Мелкая правка: 'видно.\n\n### Нако' -> 'видно.\n\n\nUPD5:\n### Нако'
en12	polosatic	2023-02-05 15:32:48	368
ru13	polosatic	2023-02-05 15:32:17	368
ru12	polosatic	2023-02-05 14:43:45	5	Мелкая правка: 'm\n{\n ld x0 = 0, x' -> 'm\n{\n int x0 = 0, x'
en11	polosatic	2023-02-05 14:43:27	5	Tiny change: 'm\n{\n ld x0 = 0, x' -> 'm\n{\n int x0 = 0, x'
en10	polosatic	2023-02-05 14:39:44	22	Tiny change: '* (b - 1)`\n' -> '* (b - 1)`, but it was obvious.\n'
ru11	polosatic	2023-02-05 14:38:59	53	Мелкая правка: 'е нашел :(\n' -> 'е нашел :(, кроме `x2 = b * (b - 1)`, что и так было очевидно.\n'
en9	polosatic	2023-02-05 14:38:03	123	Tiny change: 'nything :(\n' -> 'nything :( except that `x2 = b * (b - 1)`\n'
ru10	polosatic	2023-02-05 14:34:51	107
ru9	polosatic	2023-02-05 14:31:11	1617	Мелкая правка: ' b ^ 2`.\n*UPD4:' -> ' b ^ 2`.\n\n*UPD4:'
en8	polosatic	2023-02-05 14:29:36	1612
en7	polosatic	2023-02-05 13:13:18	14	Tiny change: 'stand why interpolation doesn't w' -> 'stand why the formula doesn't w'
ru8	polosatic	2023-02-05 13:12:54	71
en6	polosatic	2023-02-05 13:11:58	81	Tiny change: 'tion doesn`t work wit' -> 'tion doesn't work wit'
en5	polosatic	2023-02-05 13:08:43	2	Tiny change: ' use `ai >= b ^ 2`.\n' -> ' use `ai > b ^ 2`.\n'
ru7	polosatic	2023-02-05 13:08:28	1	Мелкая правка: 'вать `ai >= b ^ 2`.\n' -> 'вать `ai > b ^ 2`.\n'
ru6	polosatic	2023-02-05 12:46:30	1100	Мелкая правка: 'n\nТеперь у могу испо' -> 'n\nТеперь я могу испо'
en4	polosatic	2023-02-05 12:41:12	1081	Tiny change: '}\n }\n return ans' -> '}\n }\n return ans'
ru5	polosatic	2023-02-04 20:45:05	204
en3	polosatic	2023-02-04 20:43:23	188	Tiny change: ', a > 1.\nUPD:' -> ', a > 1.\n\nUPD:'
ru4	polosatic	2023-02-04 18:38:00	100
en2	polosatic	2023-02-04 18:37:15	106
ru3	polosatic	2023-02-04 18:36:35	4	Мелкая правка: ') = 2 * a * a - 4`, a >' -> ') = 2 * a ^ 2 - 4`, a >'
ru2	polosatic	2023-02-04 18:24:17	107	Мелкая правка: 'for the nuber of rig' -> 'for the number of rig'
en1	polosatic	2023-02-02 16:59:49	2162	Initial revision for English translation
ru1	polosatic	2023-02-02 16:29:41	2268	Первая редакция (опубликовано)

Попытки проинтерполировать

Что я хотел бы узнать:

Наконец-то нашёл формулу и решение за O(min(a, b) ^ 6)

Определения

History

Наконец-то нашёл формулу и решение за `O(min(a, b) ^ 6)`