Лучший алгоритм разложения в ряд рациональной функции
Мне нужно кодировать функцию в C ++, которая эффективно находит коэффициенты ряда Тейлора данной рациональной функции (P (x) / Q (x)).
Параметрами функции будут степень многочленов (равных по знаменателю и знаменателю), два массива с коэффициентами многочленов и числом слагаемых в разложении.
Моя идея заключалась в следующем.
Рассмотреть идентичность
P(x) / Q(x) = R(x) + ...
куда R(x) многочлен с числом слагаемых, равным количеству коэффициентов, которые мне нужно найти. Тогда я могу умножить обе стороны с Q(x) и получить
P(x) = R(x) * Q(x)
R(x) * Q(x) - P(x) = 0
Следовательно, все коэффициенты должны быть равны нулю. Это система уравнений, которая имеет O (N ^ 3) алгоритм для решения. O (N ^ 3) не так быстро, как я хотел.
Есть ли более быстрый алгоритм?
Я знаю, что коэффициенты рядов удовлетворяют линейному рекуррентному соотношению.
Это заставляет меня думать, что На) алгоритм возможен.
Решение
Алгоритм, который я собираюсь описать, математически обоснован формальный ряд власти. Каждая функция с рядом Тейлора имеет формальный степенной ряд. Обратное неверно, но если мы сделаем арифметику для функций с рядами Тейлора и получим функцию с рядом Тейлора, то мы можем сделать ту же арифметику с формальными степенными рядами и получить тот же ответ.
Алгоритм длинного деления для формальных степенных рядов подобен длинное деление алгоритм, который вы, возможно, выучили в школе. Я продемонстрирую это на примере (1 + 2 x)/(1 - x - x^2), который имеет коэффициенты, равные Номера Лукаса.
Знаменатель должен иметь ненулевой постоянный член. Мы начнем с написания числителя, который является первым остатком.
--------
1 - x - x^2 ) 1 + 2 x
Делим остаток младшего разряда (
1) постоянным членом знаменателя (1) и поставить частное сверху.
1
--------
1 - x - x^2 ) 1 + 2 x
Теперь мы умножаем 1 - x - x^2 от 1 и вычесть его из текущего остатка.
1
--------
1 - x - x^2 ) 1 + 2 x
1 - x - x^2
-------------
3 x + x^2
Сделай это снова.
1 + 3 x
--------
1 - x - x^2 ) 1 + 2 x
1 - x - x^2
---------------
3 x + x^2
3 x - 3 x^2 - 3 x^3
-------------------
4 x^2 + 3 x^3
И опять.
1 + 3 x + 4 x^2
----------------
1 - x - x^2 ) 1 + 2 x
1 - x - x^2
---------------
3 x + x^2
3 x - 3 x^2 - 3 x^3
-------------------
4 x^2 + 3 x^3
4 x^2 - 4 x^3 - 4 x^4
---------------------
7 x^3 + 4 x^4
И опять.
1 + 3 x + 4 x^2 + 7 x^3
------------------------
1 - x - x^2 ) 1 + 2 x
1 - x - x^2
---------------
3 x + x^2
3 x - 3 x^2 - 3 x^3
-------------------
4 x^2 + 3 x^3
4 x^2 - 4 x^3 - 4 x^4
---------------------
7 x^3 + 4 x^4
7 x^3 - 7 x^4 - 7 x^4
---------------------
11 x^4 + 7 x^5
Отдельные подразделения были довольно скучными, потому что я использовал делитель с ведущей 1но если бы я использовал, скажем, 2 - 2 x - 2 x^2, тогда все коэффициенты в частном будут делиться на 2,
Другие решения
Это можно сделать в O(n log n) время для произвольного P а также Q степени n, Точнее это можно сделать в M(n), где M(n) является сложностью полиномиального умножения, которое само по себе может быть сделано в O(n log n),
Во-первых, первый n члены разложения в ряд можно рассматривать просто как многочлен степени n-1,
Предположим, вы заинтересованы в первом n условия расширения серии P(x)/Q(x), Существует алгоритм, который будет вычислять обратное Q в M(n) время, как определено выше.
обратный T(x) из Q(x) удовлетворяет T(x) * Q(x) = 1 + O(x^N), То есть T(x) * Q(x) это точно 1 плюс некоторый термин ошибки, все коэффициенты которого идут после первого n Условия, которые нас интересуют, поэтому мы можем просто отбросить их.
Сейчас P(x) / Q(x) это просто P(x) * T(x) это просто еще одно полиномиальное умножение.
Вы можете найти реализацию, которая вычисляет вышеупомянутое обратное в моей библиотеке с открытым исходным кодом. Altruct. Увидеть series.h файл. Предполагая, что у вас уже есть метод, который вычисляет произведение двух полиномов, код, который вычисляет обратное, имеет длину около 10 строк (вариант «разделяй и властвуй»).
Фактический алгоритм выглядит следующим образом:
Предполагать Q(x) = 1 + a1*x + a2*x^2 + ..., Если a0 не является 1Вы можете просто разделить Q(x) а позже его обратное T(x) с a0,
Предположим, что на каждом этапе у вас есть L условия обратного, так что Q(x) * T_L(x) = 1 + x^L * E_L(x) за какую-то ошибку E_L(x), Первоначально T_1(X) = 1, Если вы включите это в выше, вы получите Q(x) * T_1(x) = Q(x) = 1 + x^1 * E_1(x) для некоторых E_1(x) что означает, что это верно для L=1, Давай теперь удвоим L на каждом шагу. Ты можешь получить E_L(x) из предыдущего шага как E_L(x) = (Q(x) * T_L(x) - 1) / x^Lили в плане реализации просто отбросьте первый L Коэффициенты произведения. Затем вы можете вычислить T_2L(x) из предыдущего шага как T_2L(x) = T_L(x) - x^L * E_L(x) * T_L(x), Ошибка будет E_2L(x) = - E_L(x)^2, Давайте теперь проверим, что шаг индукции выполнен.
Q(x) * T_2L(x)
= Q(x) * (T_L(x) - x^L * E_L(x) * T_L(x))
= Q(x) * T_L(x) * (1 - x^L * E_L(x))
= (1 + x^L * E_L(x)) * (1 - x^L * E_L(x))
= 1^2 - (x^L * E_L(x))^2
= 1 + x^2L * E_2L(x)
Что и требовалось доказать
Я почти уверен, что невозможно вычислить полиномиальное деление более эффективно, чем умножение, и, как вы можете видеть в следующей таблице, этот алгоритм только в 3 раза медленнее, чем одно умножение:
n mul inv factor
10^4 24 ms 80 ms 3,33x
10^5 318 ms 950 ms 2,99x
10^6 4.162 ms 12.258 ms 2,95x
10^7 101.119 ms 294.894 ms 2,92x
Если вы внимательно посмотрите на систему, которую вы получите с вашим планом, вы увидите, что она уже диагональна и не требует решения O (n ^ 3). Он просто вырождается в линейную рекурсию (P [], Q [] и R [] — коэффициенты соответствующих полиномов):
R[0] = P[0]/Q[0]
R[n] = (P[n] - sum{0..n-1}(R[i] * Q[n-i]))/Q[0]
Поскольку Q является полиномом, сумма имеет не более deg(Q) условия (таким образом, требуется постоянное время для расчета), делая общую сложность асимптотически линейной. Вы также можете взглянуть на матричное представление рекурсии для (возможно) лучшей асимптотики.