Всем привет. Я написал программу, которая использует эту теорему для интерполяции функции как многочлена. Сначала я расскажу об устройстве программы, а потом приведу несколько способов применения.
Код получился длинный, но пользоваться им просто.
Как использовать программу
Сначала нужно задать значения констант. 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, для того чтобы искусственно убрать старшие мономы.
Время
- Самая долгая часть предподсчета — вычисление F_CACHE — работает за O(d * O(f))
- Каждый из d запусков функции gen перебирает каждую из O(d) точек за O(N)
- Для каждого монома с ненулевым коэффициентом мы считаем его значение в каждой из 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.
