Pregunta El mejor algoritmo para la expansión en serie de la función Rational


Necesito codificar la función en C ++ que encuentre eficientemente los coeficientes de la serie de Taylor de la función racional dada (P (x) / Q (x)).

Los parámetros de función serán la potencia de los polinomios (igual en el nominador y el denominador), dos matrices con coeficientes de polinomios y número de términos en expansión.

Mi idea fue seguir. Considerar identidad

P(x) / Q(x) = R(x) + ...

Dónde R(x) es un polinomio con un número de términos igual al número de coeficientes que necesito encontrar. Entonces puedo multiplicar ambos lados con Q(x) y obten

P(x) = R(x) * Q(x)

R(x) * Q(x) - P(x) = 0

Por lo tanto, todos los coeficientes deben ser cero. Este es un sistema de ecuaciones que tienen O (n ^ 3) algoritmo para resolver. O (n ^ 3) No es tan rápido como yo quería.

¿Hay algún algoritmo más rápido?

Sé que los coeficientes de las series satisfacen la relación de recurrencia lineal. Esto me hace pensar que En) Algoritmo es posible.


5
2018-04-15 19:19


origen


Respuestas:


El algoritmo que estoy a punto de describir está justificado matemáticamente por serie de poder formal. Cada función con una serie de Taylor tiene una serie de poder formal. Lo contrario no es cierto, pero si hacemos aritmética en funciones con series de Taylor y obtenemos una función con una serie de Taylor, entonces podemos hacer la misma aritmética con series de poder formales y obtener la misma respuesta.

El algoritmo de división larga para series formales de potencia es como el división larga Algoritmo que puedes haber aprendido en la escuela. Lo demostraré en el ejemplo (1 + 2 x)/(1 - x - x^2), que tiene coeficientes iguales a los Números de Lucas.

El denominador debe tener un término constante distinto de cero. Comenzamos escribiendo el numerador, que es el primer residuo.

             --------
1 - x - x^2 ) 1 + 2 x

[ Dividimos el término de orden más bajo del residual (1) por el término constante del denominador (1) y poner el cociente arriba.

              1
             --------
1 - x - x^2 ) 1 + 2 x

Ahora nos multiplicamos 1 - x - x^2 por 1 y restarlo del residual actual.

              1
             --------
1 - x - x^2 ) 1 + 2 x
              1 -   x - x^2
              -------------
                  3 x + x^2

Hazlo otra vez.

              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

Y otra vez.

              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

Y otra vez.

              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

Las divisiones individuales eran un poco aburridas porque utilicé un divisor con un líder 1, pero si hubiera usado, digamos, 2 - 2 x - 2 x^2, entonces todos los coeficientes en el cociente se dividirían por 2.


5
2018-04-15 20:24



Esto se puede hacer en O(n log n) tiempo para arbitrario P y Q de grado n. Más precisamente, esto se puede hacer en M(n), dónde M(n) es la complejidad de la multiplicación polinomial que se puede hacer en O(n log n).

Primero de, el primero n Los términos de una expansión de la serie se pueden ver simplemente como un polinomio de grado. n-1.

Supongamos que usted está interesado en la primera n términos de la serie de expansión de P(x)/Q(x). Existe un algoritmo que computará el inverso de Q en M(n) tiempo como se define arriba.

Inverso T(x) de Q(x) satisface T(x) * Q(x) = 1 + O(x^N). Es decir. T(x) * Q(x) es precisamente 1 más un término de error cuyos coeficientes todos vienen después de la primera n términos que nos interesan, entonces podemos simplemente soltarlos.

Ahora P(x) / Q(x)es simple P(x) * T(x) que es solo otra multiplicación polinomial.

Puede encontrar una implementación que calcula el inverso mencionado anteriormente en mi biblioteca de código abierto Altructo. Ver el series.h archivo. Suponiendo que ya tiene un método que calcula el producto de dos polinomios, el código que calcula el inverso tiene una longitud de aproximadamente 10 líneas (una variante de divide y vencerás).

El algoritmo real es el siguiente: Asumir Q(x) = 1 + a1*x + a2*x^2 + .... Si a0 no es 1, simplemente puedes dividir Q(x) y luego su inverso T(x) con a0. Asume que a cada paso tienes L términos de lo inverso para que Q(x) * T_L(x) = 1 + x^L * E_L(x) por algun error E_L(x). Inicialmente T_1(X) = 1. Si enchufas esto en lo anterior obtendrás Q(x) * T_1(x) = Q(x) = 1 + x^1 * E_1(x) para algunos E_1(x) lo que significa que esto vale para L=1. Ahora doble L en cada paso. Puedes obtener E_L(x) del paso anterior como E_L(x) = (Q(x) * T_L(x) - 1) / x^L, o en la implementación, simplemente suelta el primero L Coeficientes del producto. A continuación, puede calcular T_2L(x) del paso anterior como T_2L(x) = T_L(x) - x^L * E_L(x) * T_L(x). El error sera E_2L(x) = - E_L(x)^2. Ahora comprobemos que el paso de inducción se mantiene.

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)

Q.E.D.

Estoy bastante seguro de que no es posible calcular la división polinomial más eficiente que la multiplicación, y como puede ver en la siguiente tabla, este algoritmo es solo 3 veces más lento que una sola multiplicación:

   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

2
2017-10-24 23:43



Si observa detenidamente el sistema que obtendría con su plan, puede ver que ya es diagonal y no requiere que se resuelva O (n ^ 3). Simplemente degenera en una recursión lineal (P [], Q [] y R [] son ​​los coeficientes de los polinomios correspondientes):

R[0] = P[0]/Q[0]
R[n] = (P[n] - sum{0..n-1}(R[i] * Q[n-i]))/Q[0]

Como Q es un polinomio, la suma no tiene más de deg(Q) términos (lo que toma tiempo constante para calcular), haciendo la complejidad general asintóticamente lineal. También puede mirar la representación matricial de la recursión para una (posiblemente) mejor asintótica.


1
2018-04-15 22:21