Интерполяция многочлена от нескольких переменных, или как угадать формулу

№	Пользователь	Рейтинг
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

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

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

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

Сначала нужно задать значения констант. 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.

Код

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define ld long double
#define fi first
#define se second
#define pb push_back
#define cok cout << (ok ? "YES\n" : "NO\n");
#define dbg(x) cout << (#x) << ": " << (x) << endl;
#define dbga(x,l,r) cout << (#x) << ": "; for (int ii=l;ii<r;ii++) cout << x[ii] << " "; cout << endl;
// #define int long long
#define pi pair<int, int>
const int N = 7, C = 1e7, MAX_DEG = 4, MAX_PRODUCT = 1e5;
const ld EPS = 1e-9, EPS_CHECK = 1e-9;
const string SEP = "  (", END = ")\n";
const bool APPROXIMATION = true;
array <string, N> names;
array <int, N> max_exp, powers, current_converted, cur_exp;
array<vector<ll>, N> POINTS;
ll DIV[N][MAX_DEG + 1][MAX_DEG + 1], PW[N][MAX_DEG + 1][MAX_DEG + 1];
ld SUM[MAX_PRODUCT];
ld F_CACHE[MAX_PRODUCT];
ll pow(ll a, int b)
{
	if (b == 0) return 1;
	if (b == 1) return a;
	ll s = pow(a, b / 2);
	s *= s;
	if (b & 1) s *= a;
	return s;
}
ld approximate(ld k)
{
	int k_ = k;
	int k__ = k_ + abs(k) / k;
	if (abs(k - k_) < EPS) return k_;
	else if (abs(k - k__) < EPS) return k__;
	else
	{
		int i = 1, j = 1;
		ld ka = abs(k);
		while (i < C && j < C)
		{
			ld p = ka * j;
			if (abs(p - i) < EPS) break;
			if (p < i) j++;
			else i++;
		}
		if (i >= C || j >= C) return k;
		if (k < 0) i = -i;
		return (ld)i / j;
	}	
}
void normalize(ld k)
{
    if (!APPROXIMATION)
    {
        cout << k << SEP;
        return;
    }
	int k_ = k;
	int k__ = k_ + abs(k) / k;
	if (abs(k - k_) < EPS) cout << k_ << SEP;
	else if (abs(k - k__) < EPS) cout << k__ << SEP;
	else
	{
		int i = 1, j = 1;
		ld ka = abs(k);
		while (i < C && j < C)
		{
			ld p = ka * j;
			if (abs(p - i) < EPS) break;
			if (p < i) j++;
			else i++;
		}
		if (i >= C || j >= C)
		{
			cout << k << SEP;
			return;
		}
		if (k < 0) i = -i;
		cout << i << "/" << j << SEP;
	}
}
struct monom
{
	array<int, N> exp;
	ld k;
	int deg;
	monom(array<int, N> v, ld k_)
	{
		k = k_;
		exp = v;
		deg = 0;
		for (int i=0;i<N;i++) deg += exp[i];
	}
	void display()
	{
		normalize(k);
		if (deg == 0) { cout << "1" << END; return;}
		bool go = 0;
		for (int i=0;i<N;i++)
		{
			if (go && exp[i]) cout << " * ";
			if (exp[i]) go = 1, cout << names[i] + "^" + to_string(exp[i]);
		}
		cout << END;
	}
	ld operator()(array<int, N> v)
	{
		ll res = 1;
		for (int i=0;i<N;i++) res *= PW[i][v[i]][exp[i]];
		return k * res;
	}
	ld getRandom(array<ll, N> v)
	{
		ld res = 1;
		for (int i=0;i<N;i++) res *= pow(v[i], exp[i]);
		return k * res;
	}
};
bool operator<(monom a, monom b)
{
	if (a.deg > b.deg) return 1;
	if (a.deg < b.deg) return 0;
	if (a.exp > b.exp) return 1;
	if (a.exp < b.exp) return 0;
	return a.k > b.k;
}
struct polynom
{
	vector<monom> st;
	void add(monom m)
	{
		if (abs(m.k) < EPS) return;
		st.pb(m);
	}
	void print() { if(st.size() == 0) {cout << "Polynom is 0\n"; return;} sort(st.begin(), st.end()); for (monom &m: st) m.display();}
	ld operator()(array<ll, N> v)
	{
		ld res = 0;
		for (auto &m: st) res += m.getRandom(v);
		return res;
	}
};
ld gen(int index=0, int current_hash=0)
{
	if (index == N)
	{
		ll div = 1;
		for (int i=0;i<N;i++) div *= DIV[i][current_converted[i]][cur_exp[i]];
		return (ld)(F_CACHE[current_hash] - SUM[current_hash]) / div;
	}
	ld res = 0;
	for (int i=0;i<=cur_exp[index];i++)
	{
		current_converted[index] = i;
		res += gen(index + 1, current_hash + i * powers[index]);
	}
	return res;
}
array<int, N> convert(int h)
{
	array<int, N> res;
	for (int i=0;i<N;i++) res[i] = h / powers[i], h -= res[i] * powers[i];
	return res;
}
array<ll, N> convert_points(int h)
{
	array<ll, N> res;
	for (int i=0;i<N;i++) res[i] = POINTS[i][h / powers[i]], h %= powers[i];
	return res;
}
polynom interpolate(ld f(array<ll, N>))
{
    int max_pow = -2e9, sum = 0, h_max = 0;
    set<int> remaining_points, st;
	polynom res;
    for (int x: max_exp) max_pow = max(max_pow, x), sum += x, h_max = h_max * (x + 1) + x;

    powers[N - 1] = 1;
    for (int i=N-2;i>-1;i--) powers[i] = powers[i + 1] * (max_exp[i + 1] + 1);

    for (int i=0;i<max_exp.size();i++) for (int j=0;j<=max_exp[i];j++) POINTS[i].pb(j);
	
    for (int i=0;i<N;i++) for (int j=0;j<=max_exp[i];j++) for (int u=0;u<=max_exp[i];u++) DIV[i][j][u] = (u ? DIV[i][j][u - 1] : 1) * (u == j ? 1 : (POINTS[i][j] - POINTS[i][u]));

    for (int i=0;i<N;i++) for (int j=0;j<=max_exp[i];j++) for (int u=0;u<=max_pow;u++) PW[i][j][u] = u ? PW[i][j][u - 1] * POINTS[i][j] : 1;

    for (int i=0;i<=h_max;i++) F_CACHE[i] = f(convert_points(i)), remaining_points.insert(i);
    st.insert(h_max);
	
    while (st.size())
	{
		int v = *st.rbegin();
		st.erase(v);
		remaining_points.erase(v);
		cur_exp = convert(v);
		ld k = gen();
		if (abs(k) > EPS)
		{
			if (APPROXIMATION) k = approximate(k);
			monom mn = monom(cur_exp, k);
			res.add(mn);
			for (int i: remaining_points) SUM[i] += mn(convert(i));
		}
		for (int i=0;i<N;i++) if (cur_exp[i]) st.insert(v - powers[i]);
	}
	return res;
}
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 += 13ll * i * j * u * i * i * u - 49ll * k * k * z * z * y + 90ll * c * u * k * x * x * x;
	return res;
}
void check(polynom p, ld(array<ll, N> f))
{
	mt19937 rnd(228);
	for (int i=0;i<10000;i++)
	{
		int t = clock();
		array<ll, N> ex;
		for (int j=0;j<N;j++) ex[j] = rnd() % 20 + 2;
		ld F = f(ex);
		ld P = p(ex);
		if (abs(F - P) > max(EPS_CHECK, EPS_CHECK * abs(F)))
		{
			cout << "Polynom is wrong, test " << i << endl;
			cout << F << endl << P << endl;
			for (int x: ex) cout << x << " ";
			cout << endl;
			return;
		}
		cout << "Test " << i << " has been passed, time = " << (ld)(clock() - t) / CLOCKS_PER_SEC << "s" << endl;
	}
	cout << "Polynom is OK" << endl;
}
signed main()
{
    cin.tie(0); ios_base::sync_with_stdio(0);
    cout << setprecision(20) << fixed;

    names = {"a", "b", "c", "d", "e", "f", "g"};
    max_exp = {4, 2, 3, 4, 2, 3, 3};
    
    polynom P = interpolate(f);
    P.print();
    //cout << "Checking polynom..." << endl;
    //check(P, f);
}

Стрессы

Если раскомментировать две последние строки в 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.

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