# Recurrencias Lineales

Decimos que una secuencia $s_0, s_1, s_2, \dots$ satisface una recurrencia lineal de grado $d$ si para toda $n \geq d$ se cumple que:

$$
s_n=\sum_{i=1}^{d} c_i \cdot s_{n-i}
$$

Donde $c_1, c_2, \dots, c_d$ son los coeficientes constantes que definen la recurrencia lineal y $s_0, s_1, \dots, s_{d-1}$ son los valores o términos iniciales.

:::info
:information_source: **Ejemplo**
La recurrencia lineal más conocida es la que define a la secuencia de Fibonacci: $F_n=F_{n-1}+F_{n-2}$, con los valores iniciales $F_0=0$ y $F_1=1$. Notemos que su grado es $2$, y los coeficientes son $c_1=1$ y $c_2=1$.
:::

Decimos que $s_k$ es el $k$-ésimo término de la recurrencia líneal. En programación competitiva, los problemas de recurrencias lineales normalmente involucran hallar el $k$-ésimo término.

En esta nota se explorarán los tres métodos más usados en programación competitiva para resolver recurrencias lineales:

* [DP](https://hackmd.io/@8iUL97mRS2mgfTEJDy5xmA/HyYcBzzIC) -- $O(dk)$.
* [Exponenciación de matrices](https://hackmd.io/m5cHAd4vS3Woi232-t8Y0Q?view#Exponenciaci%C3%B3n-de-Matrices) -- $O(d^3log(k))$.
* [Polinomio característico](https://hackmd.io/m5cHAd4vS3Woi232-t8Y0Q?view#Exponenciaci%C3%B3n-de-Matrices) -- $O(d^2log(k))$ o $O(dlog(d)log(k))$.

Por lo general, los problemas piden encontrar el $k$-ésimo término con valores de $k$ que pueden ir hasta $10^9$ o $10^{18}$, por lo que preferimos los métodos cuya complejidad dependen del logaritmo de $k$.

:::warning
**Nota**
Las recurrencias descritas anteriormente son llamadas **recurrencias lineales homogéneas con coeficientes constantes**, y pueden existir recurrencias lineales de otro tipo. Pero en estas notas nos centraremos en las recurrencias lineales de la forma $s_n=\sum_{i=1}^{d} c_i \cdot s_{n-i}$, que son los más comunes en problemas.
:::

## DP
La solución más simple es llevar un arreglo en el cual se van almacenando todos los términos uno por uno:

for (int n = d; n <= k; n++)
    for (int i = 1; i <= d; i++)
        v[n] += c[i] * v[n - i];

El arreglo $v$ en el código de arriba debería tener tamaño $k+1$. Sin embargo, los últimos $d$ términos para calcular $s_n$, por lo tanto podemos ahorrar memoria si sólo mantenemos los últimos $d$ elementos:

for (int n = d; n <= k; n++) {
    int sn = 0;
    for (int i = 1; i <= n; i++)
        sn += c[i] * v[d - i];
    
    // Recorrer los elementos una posición a la izquierda.
    for (int i = 0; i < d - 1; i++)
        v[i] = v[i + 1];
    v[d - 1] = sn;
}

El arreglo $v$ tiene tamaño $d$. Al principio, este arreglo contiene los términos iniciales: $v=[s_0,s_1,\ldots,s_{d-1}]$.

Justo antes de calcular el $n$-ésimo término, el arreglo contiene los términos desde el $(n-d)$-ésimo hasta el $(n-1)$-ésimo: $v=[s_{n-d}, s_{n-d+1},\ldots,s_{n-1}]$.

Al final de cada iteración, se quita el término más antiguo del arreglo, se recorre el resto a la izquierda y se almacena $s_n$ en la última posición: $v=[s_{n-d+1}, s_{n-d+2}, \ldots, s_n]$. Así se siguen mantieniendo los últimos $d$ términos en cada iteración.

Ambos códigos tienen complejidad en tiempo $O(kd)$, por lo que sólo funcionan cuando $k$ y $d$ no son muy grandes.

## Exponenciación de Matrices
Ahora vamos a ver un método que podemos utilizar cuando $k$ es grande. Este método utiliza multiplicación de matrices, por lo tanto es necesario estar familiarizado con esta operación y saber implementarla.

En el segundo código de la DP vimos que es suficiente mantener los últimos $d$ términos. Cada que calculamos un nuevo término, pasamos de tener un arreglo igual a $[s_{n-d}, s_{n-d+1},\ldots,s_{n-1}]$ a uno igual a $[s_{n-d+1}, s_{n-d+2}, \ldots, s_n]$. Si lo vemos en términos de algebra lineal, podemos decir que tenemos un vector $v$ al que le aplicamos una transormación:

$$
\begin{bmatrix}
s_{n-1}\\s_{n-2}\\ \vdots \\ s_{n-d+1} \\ s_{n-d}
\end{bmatrix}
\xrightarrow{transformación}
\begin{bmatrix}
s_{n}\\s_{n-1}\\ \vdots \\ s_{n-d+2} \\ s_{n-d+1}
\end{bmatrix}
$$

Justamente como $s_n$ es una combinación lineal de los $d$ términos anteriores, entonces se trata de una transformación lineal. Esto quiere decir que podemos utilizar una matriz (que llamaremos $T$) y multiplicarla por el vector $v$ para obtener el nuevo vector. Es decir, una matriz $T$ tal que:

$$
\begin{bmatrix}
t_{1,1}&t_{1,2}&\dots&t_{1,d-1}&t_{1,d}\\
t_{2,1}&t_{2,2}&\dots&t_{2,d-1}&t_{2,d}\\
\vdots&\vdots&\ddots&\vdots&\vdots\\
t_{d-1,1}&t_{d-1,2}&\dots&t_{d-1,d-1}&t_{d-1,d}\\
t_{d,1}&t_{d,2}&\dots&t_{d,d-1}&t_{d,d}
\end{bmatrix} \times
\begin{bmatrix}
s_{n-1}\\s_{n-2}\\ \vdots \\ s_{n-d+1} \\ s_{n-d}
\end{bmatrix}=
\begin{bmatrix}
s_{n}\\s_{n-1}\\ \vdots \\ s_{n-d+2} \\ s_{n-d+1}
\end{bmatrix}
$$

¿Cómo encontramos los valores $t_{i,j}$ de la matriz para que se realice la transformación deseada?

Observemos que, por como funciona la multiplicación de matrices, el primer elemento del vector del lado derecho será el producto punto de la primera fila de la matriz por el vector $v$, el segundo elemento será el producto punto de la segunda fila por el vector $v$, y así. Es decir, debemos elegir los valores de $t_{i,j}$ tal que:

$$
\begin{align}
t_{1,1}s_{n-1}+t_{1,2}s_{n-2}+\ldots+t_{1,d}s_{n-d} &= s_n \\
t_{2,1}s_{n-1}+t_{2,2}s_{n-2}+\ldots+t_{2,d}s_{n-d} &= s_{n-1} \\
\vdots \\
t_{d,1}s_{n-1}+t_{d,2}s_{n-2}+\ldots+t_{d,d}s_{n-d} &= s_{n-d+1}
\end{align}
$$

Para elegir los valores $t_{1,1},t_{1,2},\ldots,t_{1,d}$ de la primera fila podemos usar la definición de recurrencia lineal: $s_n=c_1s_{n-1}+c_2s_{n-2}+\ldots+c_ss_{n-s}$. Es decir, $t_{1,j}=c_j$.

Elegir los valores de las demás filas es fácil, ya que los términos $s_{n-1}, s_{n-2}, \dots, s_{n-d+1}$ ya existen en el vector del lado izquierdo. Por ejemplo, para obtener $s_{n-1}$ en la segunda fila, basta con que hagamos $1$ el elemento de la segunda fila de la matriz que se multiplica por $s_{n-1}$, y los demás los hacemos $0$, es decir, $t_{2,1}=1$ y $t_{2,j}=0$ para $j \neq 1$; en la tercera fila queremos obtener $s_{n-2}$, por lo tanto debemos hacer $t_{3,2}=1$ y $t_{3,j}=0$ para $j \neq 2$; y así. De forma general, hacemos $t_{i,i-1}=1$ a partir de la segunda fila y a todos los demás elementos los hacemos $0$.

Ahora nos queda:

$$
\begin{bmatrix}
c_1&c_2&\dots&c_{d-1}&c_{d}\\
1&0&\dots&0&0\\
0&1&\dots&0&0\\
\vdots&\vdots&\ddots&\vdots&\vdots\\
0&0&\dots&0&0\\
0&0&\dots&1&0\\
\end{bmatrix} \times
\begin{bmatrix}
s_{n-1}\\s_{n-2}\\s_{n-3}\\ \vdots \\ s_{n-d+1} \\ s_{n-d}
\end{bmatrix}=
\begin{bmatrix}
c_1s_{n-1}+c_2s_{n-2}+\ldots+c_ds_{n-d}\\
1 \cdot s_{n-1}+0\cdot s_{n-2}+\ldots+0\cdot s_{n-d}\\
0\cdot s_{n-1}+1\cdot s_{n-2}+\ldots+0\cdot s_{n-d}\\
\vdots \\
0\cdot s_{n-1}+\ldots+0\cdot s_{n-d+1}+0\cdot s_{n-d}\\
0\cdot s_{n-1}+\ldots+1\cdot s_{n-d+1}+0\cdot s_{n-d}\\
\end{bmatrix}=
\begin{bmatrix}
s_n\\s_{n-1}\\s_{n-2}\\ \vdots \\ s_{n-d+2} \\ s_{n-d+1}
\end{bmatrix}
$$

Llamaremos a la matriz $T$ que hallamos la **matriz de transición**. Sea $v_0$ el vector que contiene a los términos iniciales de la recurrencia. Veamos que si multiplicamos $T$ por $v_0$, podemos hallar $s_d$ en la primera fila del vector resultante:

$$
\begin{bmatrix}
c_1&c_2&\dots&c_{d-1}&c_{d}\\
1&0&\dots&0&0\\
0&1&\dots&0&0\\
\vdots&\vdots&\ddots&\vdots&\vdots\\
0&0&\dots&0&0\\
0&0&\dots&1&0\\
\end{bmatrix} \times
\begin{bmatrix}
s_{d-1}\\s_{d-2}\\s_{d-3}\\ \vdots \\ s_{1} \\ s_{0}
\end{bmatrix}=
\begin{bmatrix}
s_d\\s_{d-1}\\s_{d-2}\\ \vdots \\ s_{2} \\ s_{1}
\end{bmatrix}
$$

Si volvemos a multiplicar $T$ por ese resultado, podemos hallar $s_{d+1}$, y así. De forma general, si multiplicamos $x$ veces a $T$ por $v_0$, entonces podemos encontrar $s_{d-1+x}$. Por lo tanto, si queremos hallar el $k$-ésimo término de la recurrencia, debemos multiplicar $k-d+1$ veces la matriz de transición por $v_0$:

$$
\underbrace{
\begin{bmatrix}
c_1&c_2&\dots&c_{d}\\
1&0&\dots&0\\
0&1&\dots&0\\
\vdots&\vdots&\ddots&\vdots\\
0&0&\dots&0\\
0&0&\dots&0\\
\end{bmatrix} \times \dots \times
\begin{bmatrix}
c_1&c_2&\dots&c_{d}\\
1&0&\dots&0\\
0&1&\dots&0\\
\vdots&\vdots&\ddots&\vdots\\
0&0&\dots&0\\
0&0&\dots&0\\
\end{bmatrix} \times
\begin{bmatrix}
c_1&c_2&\dots&c_{d}\\
1&0&\dots&0\\
0&1&\dots&0\\
\vdots&\vdots&\ddots&\vdots\\
0&0&\dots&0\\
0&0&\dots&0\\
\end{bmatrix}
}_{k-d+1 \; veces} \times
\begin{bmatrix}
s_{d-1}\\s_{d-2}\\s_{d-3}\\ \vdots \\ s_{n-1} \\ s_{0}
\end{bmatrix}=
\begin{bmatrix}
s_k\\s_{k-1}\\s_{k-2}\\ \vdots \\ s_{k-d+2} \\ s_{k-d+1}
\end{bmatrix}
$$

La multipliación de matrices es asociativa, por lo tanto podríamos multiplicar primero la matriz $T$ por sí misma $k-d+1$ veces, es decir, elevarla a la $(k-d+1)$-ésima potencia, y después multiplicarla por $v_0$:

$$
\begin{bmatrix}
c_1&c_2&\dots&c_{d}\\
1&0&\dots&0\\
0&1&\dots&0\\
\vdots&\vdots&\ddots&\vdots\\
0&0&\dots&0\\
0&0&\dots&0\\
\end{bmatrix}^{k-d+1} \times
\begin{bmatrix}
s_{d-1}\\s_{d-2}\\s_{d-3}\\ \vdots \\ s_{n-1} \\ s_{0}
\end{bmatrix}=
\begin{bmatrix}
s_k\\s_{k-1}\\s_{k-2}\\ \vdots \\ s_{k-d+2} \\ s_{k-d+1}
\end{bmatrix}
$$

Y por lo mismo de que la multiplicación de matrices es asociativa, podemos utilizar exponenciación binaria para hallar $M^{k-d+1}$. De esta forma, sólo realizamos $O(log(k))$ multiplicaciones de matrices.

La multiplicación de dos matrices cuadradas de dimensión $d$ nos cuesta $O(d^3)$, por lo tanto la complejidad total nos queda $O(d^3log(k))$. Este método funciona bien cuando $k$ es muy grande (hasta $10^{18}$) y $d$ pequeña (apróx. $d \leq 200$). Si la $d$ es grande, necesitamos el [método que usa el polinomio característico](https://hackmd.io/m5cHAd4vS3Woi232-t8Y0Q?view#Polinomio-caracter%C3%ADstico).

### Otros tipos de recurrencias lineales
*En construcción*

## Polinomio característico
Este método hace uso de polinomios para encontrar el $n$-ésimo término de una recurrencia lineal, por lo que es preferible estar familiarizado con las [operaciones de polinomios](https://hackmd.io/@alaneos777/BJDfsAkVd) para entender la idea (aunque también es posible utilizar el [código](https://github.com/alaneos777/reference/blob/0108f4f8fe9c7f3072dbdeeded4320586d36ac78/recurrence.cpp#L12) como caja negra).

### Preliminares
Revisemos algunos conceptos antes de ver la idea detrás de este método.

#### Algoritmo de la división
El algoritmo de la división, también conocido como división euclidiana o euclídea, es un teorema que establece que, dados dos enteros $a$ y $b$ ($b \neq 0$), existen dos enteros únicos $q$ y $r$ tal que $a=bq+r$ y $0 \leq r < |b|$.

Este teorema también se puede extender a los polinomios. Dados dos polinomios $A$ y $B$ ($B \neq 0$), existen otros dos polinomios $Q$ y $R$ tal que $A=BQ+R$ y $grad(R)<grad(B)$, donde $grad(P)$ es el grado del polinomio $P$ y $deg(0)=-\infty$. Llamaremos a $Q$ el cociente y $R$ el residuo de la división euclideana entre $A$ y $B$.

:::spoiler **Demostración**
Si existe un polinomio $S$ tal que $A=SB$, entonces simplemente ponemos $Q=S$ y $R=0$. En otro caso, sea $R=A-BQ$ el menor polinomio de la forma $A-BS$ donde $S$ es un polinomio.

Demostraremos que $deg(R)<deg(B)$ por contradicción. Supongamos que $deg(R) \geq deg(B)$. Sea $r_n x^n$ el término diferente de $0$ con mayor exponente en $R$, y sea $b_m x^m$ dicho termino en $B$. Entonces $R-r_nb_m^{-1}x^{n-m}B=A-(Q+r_nb_m^{-1}x^{n-m})B$ tiene menor grado que $R$ y es de la forma $A-SB$, lo que es una contradicción.
:::

Ahora podemos definir la operación **módulo** en polinomios, la cual denotaremos con $\%$. Si $R$ es el residuo de la divisón de $A$ entre $B$, entonces $A \% B = R$.

:::info
**Ejemplo**
Consideremos $A(x)=2x^4+x^3-3x+5$ y $B(x)=x^2+1$. Tenemos que $Q(x)=2x^2+x-2$ y $R(x)=-4x+7$ cumplen que $A=BQ+R$ y $grad(R)<grad(B)$ (ya que $grad(R)=1$ y $grad(B)=2$).

Por lo tanto $(2x^4+x^3-3x+5) \% (x^2+1) = -4x+7$.
:::

El cociente y el residuo se pueden hallar en tiempo $O(n^2)$ utilizando la división larga en polinomios, o en tiempo $O(nlog(n))$ utilizando el método descrito en [estas notas de alaneos777](https://hackmd.io/@alaneos777/BJDfsAkVd#Cociente-y-residuo).

#### Polinomio característico
El polinomio característico de una recurrencia lineal de grado $d$ definida por los coeficientes $c_1, c_2, \dots c_d$ es:

$$
P(x)=x^d - \sum_{i=1}^d c_i x^{d-i}
$$

:::info
**Ejemplo**
El polinomio característico de la recurrencia lineal de grado tres $s_n=5s_{n-1}-3s_{n-2}+s_{n-3}$ es:

$$
P(x)=x^3-5x^2+3x-1
$$
:::

Este polinomio se puede usar para resolver recurrencias lineales. Por ejemplo, se podría encontrar una fórmula cerrada, como es el caso de la [Fórmula de Binet](https://artofproblemsolving.com/wiki/index.php/Binet%27s_Formula) para la secuencia de Fibonacci. También permite representar el $n$-ésimo término como una combinación lineal de los términos iniciales, que es lo mismo que ocurre en la exponenciación de matrices, pero aquí es en mejor tiempo.

#### Función $G(P)$
Definamos la función $G(P)$ asociada a una secuencia $s_0, s_1, s_2, \dots$, la cual recibe un polinomio $P(x)=p_0+p_1x+\dots+p_nx^n$, de la siguiente forma:

$$
G(P)=\sum_{i=0}^n s_ip_i
$$

Podemos ver a esta función como una especie de *producto punto* entre los coeficientes del polinomio $P$ y los términos de la secuencia.

:::info
**Ejemplo**
Consideremos la secuencia de Fibonacci $0, 1, 1, 2, 3, 5, \dots$ y $P(x)=5+x-2x^2+x^3+7x^4-8x^5$. Tenemos que

$$
G(P)=0 \cdot 5 + 1 \cdot 1 + 1 \cdot (-2) + 2 \cdot 1 + 3 \cdot 7 + 5 \cdot (-8) = -18
$$
:::

Sean $P$ y $Q$ dos polinomios, y $k$ un número real, la función $G$ tiene las siguientes propiedades:
1. $G(kP)=kG(P)$
2. $G(P+Q)=G(P) + G(Q)$

:::spoiler **Demostración**
1. $G(kP)=\sum_{i} s_i (kp_i)=k\sum_{i} s_i p_i=kG(P)$
2. $G(P+Q)=\sum_{i} s_i (p_i+q_i)=\sum_{i} s_i p_i+\sum_{i} s_i q_i=G(P)+G(Q)$
:::

Ahora supongamos que la secuencia $s_0, s_1, s_2, \dots$ satisface una recurrencia lineal definida por los coeficientes $c_1, c_2, \dots, c_d$. Si $P$ es el polinomio caracteristico de dicha recurrencia lineal, veamos cuánto vale $G(P)$.

$$
\begin{align}
G(P) &= G(x^d-\sum_{i=1}^dc_ix^{d-i}) \\
&= G(x^d-c_1x^{d-1}-c_2x^{d-2}-\ldots-c_d) \\
&= s_d-s_{d-1}c_1-s_{d-2}c_2-\ldots - s_0c_d \\
&= s_d-\sum_{i=1}^d s_{d-i}c_i \\
&= s_d-s_d \\
&= 0
\end{align}
$$

Así podemos observar que $G(P)=0$ para el polinomio característico de una recurrencia lineal. También se cumplirá que $G(x^kP)=0$ para cualquier $k\geq0$.

:::spoiler **Demostración**
$$
\begin{align}
G(x^kP) &= G(x^{d+k}-\sum_{i=1}^dc_ix^{d-i+k}) \\
&= G(x^{d+k}-c_1x^{d-1+k}-c_2x^{d-2+k}-\ldots-c_dx^k) \\
&= s_{d+k}-s_{d-1+k}c_1-s_{d-2+k}c_2-\ldots - s_kc_d \\
&= s_{d+k}-\sum_{i=1}^d s_{d+k-i}c_i \\
&= s_{d+k}-s_{d+k} \\
&= 0
\end{align}
$$
:::

La última propiedad que veremos sobre $G$ es que si $P$ es el polinomio caracteristico y $Q$ es cualquier polinomio, entonces $G(QP)=0$

:::spoiler **Demostración**
$$
\begin{align}
G(QP) &= G((q_0+q_1x+\ldots+q_nx^n)P) \\
&= G((q_0P)+(q_1xP)+\ldots+(q_nx^nP)) \\
&= G(q_0P)+G(q_1xP)+\ldots+G(q_nx^nP) \\
&= q_0G(P)+q_1G(xP)+\ldots+q_nG(x^nP) \\
&= q_0\cdot 0+q_1\cdot 0+\ldots+q_n\cdot 0 \\
&= 0
\end{align}
$$
:::

### Método para encontrar el $k$-ésimo término de una recurrencia lineal
Ya tenemos todo lo necesario para encontrar el $k$-ésimo término de una recurrencia lineal de grado $d$ utilizando el polinomio característico.

Queremos evaluar $G(x^k)$, ya que $G(x^k)=s_k \cdot 1=s_k$. Puede parecer que $G$ no nos ayuda mucho para hallar $s_k$, precisamente porque necesitamos a $s_k$ para evaluar $G(x^k)$. ¿De qué forma podemos evaluar $G(x^k)$ sin tener que calcular $s_k$ antes?

El algoritmo de la división nos ayuda a encontrar un polinomio $R$ de grado menor a $d$ a partir de $x^k$:

$$
x^k=Q(X)\cdot P(X) + R(x)
$$

Ya vimos que $G(QP)=0$ si $P$ es el polinomio característicon de la recurrencia. Esa propiedad nos ayuda a quedarnos solo con $R(x)$:

$$
\begin{align}
s_k &= G(x^k) \\
&= G(Q(x)P(x)+R(x)) \\
&= G(Q(x)P(x)) + G(R(x)) \\
&= G(R(x))
\end{align}
$$

Así pues, vemos que $G(R(x))=s_k$, donde $R(x)=x^k \% P(x)$. Como $P(x)$ tiene grado $d$, entonces $R(d)$ tiene grado menor a $d$, por lo que para calcular $G(R(x))$ necesitamos únicamente los términos iniciales de la secuencia.

Si ya tenemos $R$, podemos calcular $s_k$ en complejidad $O(d)$, así que el cuello de botella en la complejidad ahora está en calcular $x^k \% P(x)$. Como se mencionó en la [sección del algoritmo de la división](https://hackmd.io/m5cHAd4vS3Woi232-t8Y0Q?view#Algoritmo-de-la-divisi%C3%B3n), la operación módulo en polinomios se puede hacer en $O(k^2)$ o $O(klog(k))$, lo que todavía no es óptimo en tiempo, ya que $k$ puede ir hasta $10^{18}$ en problemas de recurrencia lineal.

Utilicemos una idea que usamos en números enteros cuando tenemos que encontrar una potencia grande módulo un entero: exponenciación binaria, pero ahora con polinomios. El código quedaría de la siguiente forma:

int resuelveRecurrenciaLineal (vector<int> c, vector<int> s, int d, int k) {
    polinomio R = {1}, U = {0, 1}; // R = 1, U = x
    
    polinomio P(d + 1);
    P[d] = 1;
    for (int i = 1; i <= d; i++)
        P[d - i] = -c[i];
    
    while (k) {
        if (k & 1)
            R = (R * U) % P;
        U = (U * U) % P;
        k >>= 1;
    }
    
    int res = 0;
    for (int i = 0; i < R.degree(); i++)
        res += R[i] * s[i];
    return res;
}

Al modular por $P$ en cada iteración de la exponenciación binaria, $R$ y $U$ tienen grado menor a $d$, por lo tanto ahora la complejidad de hallar $x^k \% P(x)$ dependen de $d$ y $log(K)$:

* Si se usa la múltiplicación bruta de polinomios y el algoritmo largo de la división para los módulos, la complejidad es $O(d^2log(k))$.
* Si se usa la múltiplicación con FFT y el algoritmo de modulación expuesto en [estas notas](https://hackmd.io/@alaneos777/BJDfsAkVd#Cociente-y-residuo), la complejidad es $O(dlog(d)log(k))$.