Блог пользователя polosatic

Автор polosatic, история, 3 недели назад, перевод, По-русски

Всем привет! Сегодня я хочу познакомить вас с новым алгоритмом: тернарным поиском на неунимодальных функциях. Я придумал ее во время ROI 2024, и благодаря этому стал призером.

Мотивация

В олпроге есть много примеров, когда нам нужно максимизировать или минимизировать значение некоторой функции $$$f$$$. Если функция унимодальна, то нет проблем с использованием тернарного поиска. Но если это неправда, то делать нечего. Конечно, если $$$f$$$ случайна, нужно проверить все ее значения. К счастью, интуиция подсказывает нам, что если функция близка к унимодальной, то мы можем использовать тернарный поиск и на большей части ее области определения.

Интуиция

Рассмотрим унимодальную функцию $$$0.03x^2$$$, добавив к ней немного шума ($$$sin (x)$$$):

Graph

Интуиция подсказывает, что если шум не очень сильный, мы все равно можем использовать троичный поиск вдали от глобального минимума.

В данном примере мы можем показать это более формально: производная функции равна $$$0.06x + cos(x)$$$, и при $$$|x| > \frac{1}{0.06} \approx 16.67$$$, добавление косинуса не меняет знак производной.

Первая оптимизация

Итак, первая идея — запустить тернарный поиск, используя не while (r - l > eps), а while (r - l > C), и затем перебрать все значения между $$$l$$$ и $$$r$$$ с некоторой точностью. В случае, когда $$$f$$$ принимает целое число, проблема точности вообще отсутствует.

Вторая оптимизация

Она упоминается в этом блоге. В нем рассказывается аналогичная идея разбиения всех значений аргументов на блоки и применения к ним троичного поиска.

Это единственная связанная с блогом идея, которую я нашел в Интернете. Я пробовал гуглить и спрашивать людей, занимающихся олпрогой, но никто из них не знал об этом раньше.

Тесты

Функция из примера скучна, поэтому рассмотрим более интересную: функция Вейерштрасса.

Приближая, я выяснил, что максимум составляет около $$$1,162791$$$.

Spoiler

Будем искать максимум на интервале (-1, 1).

Code

Это дает нам $$$1,12881$$$. Изменение $$$eps$$$ лишь немного меняет это значение.

Давайте разделим аргументы на блоки. Поскольку аргументы вещественные, мы не будем явно разбивать их, а возьмем максимум в некотором диапазоне вблизи аргумента.

Blocks

Это дает $$$1,15616$$$, что весьма неплохо. Мы можем еще улучшить ответ, взяв максимум из всех значений $$$f$$$, которые мы когда-либо вычисляли:

Take maximum

Это дает нам $$$1,16278$$$, что очень близко к ожидаемым $$$1,162791$$$. Кажется, что мы нашли хороший метод. Но давайте копнём глубже.

Третья оптимизация

Давайте изменим константу $$$3$$$ в коде. Назовем ее $$$C$$$. Опытным олпрогерам известно, что часто хорошо выбрать $$$C$$$ равной $$$2$$$ (бинарный поиск по производной) или $$$\frac{\sqrt5+1}{2}$$$ (терпоиск с золотым сечением). Поскольку на каждой итерации мы отбрасываем $$$\frac{1}{C}$$$ часть нашего интервала, вероятность того, что максимум окажется внутри отрезанной части, уменьшается с увеличением C.

Выберем $$$C = 100$$$ и запустим тернарный поиск. Как мы уже знаем, взятие максимума из всех вызовов функций — очень хорошая оптимизация, поэтому я сразу ее добавлю.

Code

Если мы запустим его при $$$C = 3$$$, то получим $$$1,13140$$$ (алгоритм тот же, что и в классическом троичном поиске, но мы берем максимум, поэтому ответ намного лучше). Теперь давайте увеличим $$$C$$$ и посмотрим, как увеличится ответ:

Мы запускаем его с $$$C = 30$$$ и получаем $$$1,16275$$$. Запускаем его с $$$C = 100$$$ и получаем... $$$1,15444$$$.

Дело в том, что увеличение $$$C$$$ не гарантирует увеличения ответа.

Четвертая оптимизация

Давайте переберем все значения $$$C$$$ от 3 до 100 и возьмем лучший ответ:

for (int C = 3; C < 100; C++) res = max(res, find_maximum(Weierstrass, C));

Это дает нам $$$1,16279$$$ и работает быстрее, чем разбиение на блоки. Чтобы сравнивать эти подходы дальше, нам нужно изменить функцию, потому что оба метода возвращают значения, которые очень близки к ответу.

Используем $$$a = 0,88$$$ и $$$b = 51$$$ для функции Вейерштрасса. Обратите внимание, что в Desmos невозможно увидеть фактический максимум функции.

Я сравню эти два подхода в Запуске Codeforces.

C < 100
C < 200
C < 400

Я попробовал объединить эти методы вместе, но результат оказался хуже, чем у обоих методов (потому что я использовал $$$1000$$$ итераций вместо $$$10000$$$ и перебор $$$C < 10$$$, чтобы влезть в ТЛ).

Вывод

На рассмотренной мной функции оба подхода достаточно хороши. На реальных соревнованиях я использовал только свой подход (четвертую оптимизацию).

Из недавних задач, 2035E - Монстр, я успешно использовал эту технику. Вот реализация для целочисленных аргументов: 288346163.

Если вы знаете, как найти максимум лучше (не методом имитации отжига), пожалуйста, напишите.

Полный текст и комментарии »

  • Проголосовать: нравится
  • +86
  • Проголосовать: не нравится

Автор polosatic, история, 3 месяца назад, По-русски

Недавно я решал эту задачу. Очевидно, тут нужны компоненты сильной связности. Я их написал и получил TL14: 277474774.

Я знал, что пушбэки в 1e6 векторов очень медленные из-за динамической реаллокации. Я также знал единственный способ исправить это: предподсчитать размер каждого вектора, и интерепретировать vector<vector> как один статический массив. Реализовывать это было неприятно, но я справился и получил TL38: 277476691

Я пытался добавить разные оптимизации: убирал последний вектор Эйлерова обхода, добавлял прагмы, но решение не проходило 38 тест. Поэтому, у меня возникло 2 вопроса для самых опытных оптимизаторов:

  1. Есть ли другой, более простой способ ускорить пушбэки в векторы, без предподсчета размера каждого вектора?

  2. Почему мое решение падает по TL и как его ускорить? Ограничения до $$$10^6$$$ должны влезать в 2.5 секунды.

UPD: #define vector basic_string — очень хорошая оптимизация. EDIT: не используйте этот define если пользуетесь vector<vector<int>> vec вместо vector<int> vec[N].

UPD2: нашел этот блог. Там объяснены нюансы.

Полный текст и комментарии »

  • Проголосовать: нравится
  • +46
  • Проголосовать: не нравится

Автор polosatic, история, 6 месяцев назад, перевод, По-русски

Недавно я решал эту задачу. Я свёл ее к более простой задаче: даны точки A, B, C, и надо узнать точку пересечения 2 прямых: AB и OC, где точка O — начало координат.

Код, 90/92 тестов пройдено

Сначала я написал пересечение двух прямых формулой (функция find()), и получил ВА на 2 тестах. Потом я написал это при помощи бинпоиска (функция findgood()) и получил AC. Далее я снова заслал решение с find(), добавив #define ld __float128, и оно зашло.

Кто-нибудь может объяснить, что за чертовщина тут произошла? Как формульное решение может быть менее точным, чем бинпоиск?

Полный текст и комментарии »

  • Проголосовать: нравится
  • +16
  • Проголосовать: не нравится

Автор polosatic, история, 8 месяцев назад, По-русски

В физике есть много читерских утверждений, вроде законов сохранения, которые помогают "доказывать" разные математические факты. Блог создан для того, чтобы люди делились известными им доказательствами. Вот что известно мне:

Нормали тетраэра

Условие
Решение

Теорема Пифагора

Условие
Решение

Точка Торричелли

Условие
Решение

Полный текст и комментарии »

  • Проголосовать: нравится
  • +15
  • Проголосовать: не нравится

Автор polosatic, история, 11 месяцев назад, По-английски

(Unfortunately, Codeforces does not support emoji, so it is here)

Полный текст и комментарии »

  • Проголосовать: нравится
  • +45
  • Проголосовать: не нравится

Автор polosatic, история, 18 месяцев назад, По-русски

Всем привет. Я написал программу, которая использует эту теорему для интерполяции функции как многочлена. Сначала я расскажу об устройстве программы, а потом приведу несколько способов применения.

Код получился длинный, но пользоваться им просто.

Как использовать программу

Сначала нужно задать значения констант. N — это количество переменных в Вашем многочлене, MAX_DEG — это максимально возможная степень переменной, в которой она может входить в какой-либо из одночленов. В функции main нужно заполнить два массива N элементами: names содержит имена всех переменных, max_exp на i-той позиции содержит максимальный показатель степени (или оценку сверху на него), который может иметь соответствующая переменная.

Обозначим d = (max_exp[0] + 1) * (max_exp[1] + 1) * ... * (max_exp[N - 1] + 1). Должно выполняться, что константа MAX_PRODUCT больше, чем d. Дальше нужно написать функцию f, которая на вход принимает array<ll, N>, а возвращает ll или ld. В моём примере, результат работы функции — целое число, но функция возвращает ld для того, чтобы избежать переполнений ll.

Код

Стрессы

Если раскомментировать две последние строки в main, то программа сама проверит получившийся многочлен на случайных тестах. Генерацию тестов нужно изменять под конкретную функцию f, иначе она может долго вычисляться на больших тестах.

Приближения

Функция из примера (и все подобные функции с N циклами) является многочленом с рациональными коэффициентами (иначе целое число на выходе мы не получим). Поэтому, в случае APPROXIMATION = true, все коэффициенты приближаются к рациональным с абсолютной погрешностью EPS при помощи функций normalize и approximate. Приближения к рациональным дробям выполняются, вероятно, не самым эффективным алгоритмом за O(числитель + знаменатель), но при небольшом количестве мономов в многочлене это недолго.

Функция стресс-тестирования считает результат вычисления многочлена корректным, если его абсолютная или относительная погрешность не больше, чем EPS_CHECK.

Как и за сколько времени это работает

Мономы мы представляем в виде массива показателей степеней переменных, которые мы хэшируем. Массив PW — предпосчёт степеней, в которые возводим числа в массиве POINTS — собственно, точки, по которым мы интерполируем. Если Вы хотите задать свои точки для интерполяции, то нужно изменить массив POINTS. Если там будут дробные числа, то в начале программы нужно заменить #define ll long long на #define ll long double. Массив DIV служит для быстрого вычисления знаменателей в формуле коэффициента.

convert(h) — получить индексы координат точки в массиве POINTS, соответствующей моному с хэшом h convert_points(h) — получить координаты точки, соответствующей моному с хэшом h.

Далее мы предподсчитываем значения функции f во всех наших точках и записываем их в массив F_CACHE. Потом мы запускаем bfs по мономам, где мы при переходе от одного монома к другому уменьшаем показатель степени одной из переменных на 1. Приходя в bfs'е к моному, мы находим коэффициент при нём при помощи функции gen. Если коэффициент ненулевой, то мы должны изменить наш многочлен для всех ещё не пройденных мономов. (Здесь мы не разделяем понятия монома и точки, так как из показателей степеней монома мы можем получить N координат точки при помощи функции convert_points(h), где h — хэш монома). Это нужно для того, чтобы выполнялось одно из условий теоремы: в многочлене не должно быть мономов старше нашего. Мы для каждой точки добавляем в массив SUM значение в этом мономе, чтобы потом в функции gen его вычесть из результата работы функции f, для того чтобы искусственно убрать старшие мономы.

Время

  1. Самая долгая часть предподсчета — вычисление F_CACHE — работает за O(d * O(f))
  2. Каждый из d запусков функции gen перебирает каждую из O(d) точек за O(N)
  3. Для каждого монома с ненулевым коэффициентом мы считаем его значение в каждой из O(d) точек за O(N)

Получили O(d * O(f) + d^2 * N + d * O(res)), где O(res) — время для вычисления полученного в результате многочлена.

Попытка оптимизировать

Скорее всего, больше всего времени будет занимать рекурсия. Её можно развернуть в цикл со стеком. Это скучно, и я решил узнать, что будет, если её развернуть просто в цикл. Давайте вместо запуска рекурсии пробежимся по всем хэшам мономов, меньших нашего. Для каждого монома проверим, является ли он младше нашего (все соответствующие показатели степеней небольше). Если младше, то добавляем к текущему коэффициенту значение дроби для этой точке. Код будет какой-то такой:

// Вместо ld k = gen();
ld k = 0;
for (int h=0;h<=v;h++)
{
    array<int, N> cur = convert(h);
    bool ok = 1;
    for (int i=0;i<N;i++) if (cur[i] > cur_exp[i]) ok = 0;
    if (ok)
    {
	ll div = 1;
        for (int i=0;i<N;i++) div *= DIV[i][cur[i]][cur_exp[i]];
        k += (ld)(F_CACHE[h] - SUM[h]) / div;
    }
}

Будет ли это быстрее? Новая реализация перебирает по 1 разу каждую пару хэшей, поэтому она работает за O(d^2 * N), как и функция gen. Теперь оценим константу. Пар хэшей существует d * (d + 1) / 2. Константа 1 / 2. Чему равна константа количества рассмотренных точек функции gen? По сути, это количество можно посчитать при помощи функции:

ld f(array<ll, N> v)
{
	auto [a, b, c, d, e, f, g] = v;
	ld res = 0;
	for (int i=0;i<a;i++)
		for (int j=0;j<b;j++)
			for (int u=0;u<c;u++)
				for (int x=0;x<d;x++)
					for (int y=0;y<e;y++)
						for (int z=0;z<f;z++)
							for (int k=0;k<g;k++)
								res += (i + 1) * (j + 1) * (u + 1) * (x + 1) * (y + 1) * (z + 1) * (k + 1);
	return res;
}

Коэффициент при a^2 * b^2 * c^2 * d^2 * e^2 * f^2 и будет нашей константой. Для нахождения этого коэффициента я воспользовался своей программой. Он оказался равен 1/128. Вообще, для N переменных он равен 1 / 2^N. То есть способ оптимизации эффективен для очень маленьких N.

Заключение

Возможно, кому-то эта программа поможет узнать формулу для какой-то функции. Также она может раскрывать скобки, что необходимо при счёте геометрии в комплексных числах. Если Вы придумали другие способы использования, то я буду рад, если Вы ими поделитесь.

При N = 1 эта программа — просто интерполяция по Лагранжу, для которой существует реализация быстрее, чем за квадрат. Возможно, кто-нибудь сможет придумать ускорение и при N > 1.

Полный текст и комментарии »

  • Проголосовать: нравится
  • +76
  • Проголосовать: не нравится

Автор polosatic, история, 21 месяц назад, По-русски

Всем привет. На задачу 1771F - Хоссам и минимум на отрезке я сделал две идентичные посылки:

TL: 193157567 AC: 193157512

Можно убедиться, что они отличаются только в одной строке с map. Не могу понять, почему так произошло, что первое решение не зашло.

Полный текст и комментарии »

  • Проголосовать: нравится
  • +8
  • Проголосовать: не нравится

Автор polosatic, история, 22 месяца назад, По-русски

Всем привет. Недавно мне в голову пришла следующая задача: Рассмотрим множество точек (x, y) с целыми координатами таких, что 0 <= x < a и 0 <= y < b. Требуется найти количество остроугольных треугольников с вершинами в этих точках.

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

Ясно, что можно написать функцию f(a, b), которая будет искать ответ и работать при этом за (ab) ^ 3. Я предположил, что она ведет себя как многочлен от двух переменных степени не более 6. Я попытался её проинтерполировать, используя эту теорему. Но у меня ничего не получилось, так как при мономах степени больше 6 интерполяция давала ненулевой коэффициент. Не получилось также с тупоугольными и прямоугольными треугольниками.

code (для прямоугольных треугольников)

Данный код узнаёт коэффициент при мономе 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

Теперь я могу использовать большие значения a и b для интерполяции.

Но всё равно кажется странным, что количество прямоугольных треугольников пропорционально (ab) ^ 2, а не (ab) ^ 3. Сейчас попробую понять, почему формула не работает для ai <= b ^ 2. EDIT: ai <= (b - 1) ^ 2

UPD4: Код, который находит формулу для f(a, b) при a > b ^ 2 и работает за O(b ^ 6):

Спойлер

Пытаясь найти закономерность в коэффициентах многочлена 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 типа:

  1. треугольники, у которых нет вершин на прямой x = a - 1

  2. треугольники, у которых есть вершины на прямой x = a - 1 и две стороны параллельны осям координат

  3. остальные треугольники

Обозначим количество треугольников третьего типа за 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 — константа.

Доказательство

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).

Преобразования
Реализация

К сожалению, c для разных b — разные, и я не смог найти между ними закономерность. Не помогли ни интерполяция, ни OEIS. Осталось несколько вещей, которые надо сделать:

  1. Выразить c через b, ведь c зависит только от b
  2. Решить задачу для a <= (b - 1) ^ 2 быстрее
  3. Решить задачу для остроугольных треугольников

Будет жёстко.

UPD6: Можно ускорить решение до O(b^5), используя эту идею

Полный текст и комментарии »

  • Проголосовать: нравится
  • +88
  • Проголосовать: не нравится