В физике есть много читерских утверждений, вроде законов сохранения, которые помогают "доказывать" разные математические факты. Блог создан для того, чтобы люди делились известными им доказательствами. Вот что известно мне:
Нормали тетраэра
Теорема Пифагора
Точка Торричелли
Всем привет. Я написал программу, которая использует эту теорему для интерполяции функции как многочлена. Сначала я расскажу об устройстве программы, а потом приведу несколько способов применения.
Код получился длинный, но пользоваться им просто.
Сначала нужно задать значения констант. 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, для того чтобы искусственно убрать старшие мономы.
Получили 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.
Всем привет. На задачу 1771F - Hossam and Range Minimum Query я сделал две идентичные посылки:
Можно убедиться, что они отличаются только в одной строке с map. Не могу понять, почему так произошло, что первое решение не зашло.
Всем привет. Недавно мне в голову пришла следующая задача: Рассмотрим множество точек (x, y) с целыми координатами таких, что 0 <= x < a и 0 <= y < b. Требуется найти количество остроугольных треугольников с вершинами в этих точках.
Ясно, что можно написать функцию f(a, b), которая будет искать ответ и работать при этом за (ab) ^ 3. Я предположил, что она ведет себя как многочлен от двух переменных степени не более 6. Я попытался её проинтерполировать, используя эту теорему. Но у меня ничего не получилось, так как при мономах степени больше 6 интерполяция давала ненулевой коэффициент. Не получилось также с тупоугольными и прямоугольными треугольниками.
Данный код узнаёт коэффициент при мономе a ^ (C1 - 1) * b ^ (C2 - 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).
Теперь я могу использовать большие значения 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 типа:
треугольники, у которых нет вершин на прямой 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 — константа.
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. Осталось несколько вещей, которые надо сделать:
c через b, ведь c зависит только от ba <= (b - 1) ^ 2 быстрееБудет жёстко.
UPD6: Можно ускорить решение до O(b^5), используя эту идею
| № | Пользователь | Рейтинг |
|---|---|---|
| 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 |
| Страны | Города | Организации | Всё → |
| № | Пользователь | Вклад |
|---|---|---|
| 1 | maomao90 | 174 |
| 2 | awoo | 164 |
| 3 | adamant | 163 |
| 4 | TheScrasse | 159 |
| 5 | nor | 157 |
| 6 | maroonrk | 156 |
| 7 | -is-this-fft- | 152 |
| 8 | Petr | 146 |
| 8 | orz | 146 |
| 10 | BledDest | 145 |
