FFT

La Transformada Rápida de Fourier (FFT) se basa en números complejos, aunque no es esencial comprenderlos a fondo para su implementación práctica.

Ahora ya lo sabemos, ¿verdad?.

Reflexión

Consideremos cómo multiplicar dos polinomios \(f,g\). Un enfoque directo es con complejidad \(\mathcal O(n\times m)\), pero no parece optimizable de manera evidente.

Podemos pensar en trabajar desde el resultado: ya que el producto es un polinomio de grado \(N=n+m\), ¿por qué no interpolarlo directamente? Esto parece correcto, pero resulta en \(\mathcal O(N^2)\), ¡incluso más lento!

Propiedad Clave

Según la interpolación de Lagrange, podemos elegir cualquier \(N\) valores para evaluar. Seleccionar \(1\sim N\) no ayuda. Aquí se requiere imaginación: usar las raíces de la unidad de orden \(N\), específicamente sus potencias.

Método

Evaluación Polinomial

Esencialmente, es multiplicar un vector por una matriz:

[\begin{bmatrix} x_{0} \ x_{1} \ x_{2} \ \vdots \ x_{N-1} \end{bmatrix} \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \ 1 & \alpha & \alpha^{2} & \cdots & \alpha^{N-1} \ 1 & \alpha^{2} & \alpha^{4} & \cdots & \alpha^{2(N-1)} \ \vdots & \vdots & \vdots & \ddots & \vdots \ 1 & \alpha^{N-1} & \alpha^{2(N-1)} & \cdots & \alpha^{(N-1)(N-1)} \end{bmatrix} = \begin{bmatrix} X_{0} \ X_{1} \ X_{2} \ \vdots \ X_{N-1} \end{bmatrix} ]Sea \(\alpha\) una raíz de la unidad de orden \(N\), entonces la matriz tiene solo \(N\) elementos distintos, pero aún es costoso calcular.

Aprovechemos una propiedad: si \(N\) es par, \(\alpha^k=-\alpha^{k+\frac{N}{2}}\). Esto permite calcular \(f(\alpha^k),f(\alpha^{k+\frac{N}{2}})\) simultáneamente, reduciendo el trabajo a \(\frac{N}{2}\) evaluaciones. Sin embargo, si evaluamos \(f(\alpha^k)\) directamente, no hay mejora.

Separemos los términos pares e impares de \(f\):

[f_0(x)=a_0+a_2x^2...+a_{2d}x^{2d} ][f_1(x)=a_1x+a_3x^3...+a_{2k+1}x^{2k+1}=x(a_1+a_3x^2...+a^{2d+1}x^{2d}) ]Al sustituir \(\alpha^k \):

[f_0(\alpha^k)=a_0+a_2\alpha^{2k}+...a_{2d}\alpha^{2kd}=f_0(\alpha^{2k}) ]Análogamente,

[f_1(\alpha^k)=\alpha^{k}f_1(\alpha^{2k}) ]Nótese que \(\alpha_{N}^{2k}\) es equivalente a \(\alpha_{\frac{N}{2}}^k\) bajo la raíz de la unidad de orden \(\frac{N}{2}\), convirtiéndose en un subproblema.

[f(\alpha^k)=f_0(\alpha^{2k})+\alpha^kf_1(\alpha^{2k})=f'(\alpha_{\frac{N}{2}}^{k})+\alpha^kf''(\alpha_{\frac{N}{2}}^{k}) ]Las funciones \(f'(x)\) y \(f''(x)\) tienen la misma forma que \(f(x)\), pero con grado reducido a la mitad, permitiendo recursión.

Además,

[f(\alpha^{k+\frac{N}{2}})=f_0(\alpha^{2k})+\alpha^{k+\frac{N}{2}}f_1(\alpha^{2k})=f'(\alpha_{\frac{N}{2}}^{k})-\alpha^kf''(\alpha_{\frac{N}{2}}^{k}) ]Así, ambos valores se calculan juntos.

Recuperación del Polinomio

Actualmente solo evaluamos en raíces de la unidad; interpolar para obtener coeficientes parece \(\mathcal O(N^2)\). Dado que la evaluación multiplica por una matriz, la recuperación usa la matriz inversa.

[\begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \ 1 & \alpha & \alpha^{2} & \cdots & \alpha^{N-1} \ 1 & \alpha^{2} & \alpha^{4} & \cdots & \alpha^{2(N-1)} \ \vdots & \vdots & \vdots & \ddots & \vdots \ 1 & \alpha^{N-1} & \alpha^{2(N-1)} & \cdots & \alpha^{(N-1)(N-1)} \end{bmatrix} ]Su inversa es especial:

[\frac{1}{N} \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \ 1 & \alpha^{N-1} & \alpha^{2(N-1)} & \cdots & \alpha^{(N-1)(N-1)} \ 1 & \alpha^{(N-2)} & \alpha^{2(N-2)} & \cdots & \alpha^{(N-2)(N-1)} \ \vdots & \vdots & \vdots & \ddots & \vdots \ 1 & \alpha^{} & \alpha^{2} & \cdots & \alpha^{(N-1)} \end{bmatrix} ]Cada \(\alpha^k\) se transforma en su conjugado \(\alpha^{N-k}\). Entonces, se aplica la evaluación inversa.

Complejidad: \(\mathcal O(n\times \log n)\).

Implementación

Debido al uso de números de punto flotante, la versión recursiva es lenta. Analizando la esencia, se puede derivar una versión iterativa.

Primero, la FFT recursiva separa términos pares e impares en cada nivel. Al llegar al caso base, es asignar valores según la inversión binaria. Esto se puede implementar con un bucle.

Luego, la fusión usa la "optimización de mariposa", calculando in-situ \(\alpha^k\) y \(\alpha^{k+\frac{N}{2}}\) para mejor comprensión.

Código de ejemplo:

#include<bits/stdc++.h>
using namespace std;
const int maxn=5e6+5;
const double dos_pi=2.0*acos(-1);
int n,m;
struct Complex{
	double real, imag;
	Complex operator + (const Complex &otro) const {
		return {real + otro.real, imag + otro.imag};
	}
	Complex operator - (const Complex &otro) const {
		return {real - otro.real, imag - otro.imag};
	}
	Complex operator * (const Complex &otro) const {
		return {real * otro.real - imag * otro.imag, real * otro.imag + imag * otro.real};
	}
} roots[maxn], polyA[maxn], polyB[maxn];

void fourierTransform(Complex arr[], int len, int dir){
	for(int i=0, j=0; i<len; i++){
		if(i < j) swap(arr[i], arr[j]);
		for(int k=len/2; j^=k, j<k; k>>=1);
	}
	roots[0] = {1,0};
	for(int step=1; step<len; step<<=1){
		double angle = 1.0/(step*2);
		Complex base = {cos(dos_pi*angle), dir*sin(dos_pi*angle)};
		for(int j=step*2-2; j>=0; j-=2){
			roots[j] = roots[j>>1];
			roots[j+1] = roots[j] * base;
		}
		for(int i=0; i<len; i+=(step<<1)){
			for(int j=i; j<i+step; j++){
				Complex temp1 = arr[j], temp2 = arr[j+step] * roots[j-i];
				arr[j] = temp1 + temp2, arr[j+step] = temp1 - temp2;
			}
		}
	}
	if(dir == -1){
		double invLen = 1.0/len;
		for(int i=0; i<len; i++){
			arr[i].real *= invLen, arr[i].imag *= invLen;
		}
	}
}
int main(){
	scanf("%d%d", &n, &m);
	for(int i=0; i<=n; i++){int val; scanf("%d", &val); polyA[i].real=val;}
	for(int i=0; i<=m; i++){int val; scanf("%d", &val); polyB[i].real=val;}
	int totalLen = n+m+1;
	if((1<<__lg(totalLen)) != totalLen) totalLen = (1<<(__lg(totalLen)+1));
	fourierTransform(polyA, totalLen, 1);
	fourierTransform(polyB, totalLen, 1);
	for(int i=0; i<totalLen; i++) polyA[i] = polyA[i] * polyB[i];
	fourierTransform(polyA, totalLen, -1);
	for(int i=0; i<=n+m; i++) printf("%d ", int(polyA[i].real+0.5));
	return 0;
}

NTT

Aplicación

Se utiliza para multiplicar polinomios con coeficientes módulo un primo.

Idea

Notamos que la raíz primitiva es análoga a la raíz de la unidad. Definimos la raíz de la unidad de orden \(N\) módulo \(m\) como \(\omega^N\), representada por \(g^{\frac{m-1}{N}}\), donde \(g\) es una raíz primitiva. Las propiedades son similares.

Dado que \(\varphi(m)=m-1\) se divide por \(2\) en cada paso, los módulos NTT suelen tener \(m-1\) con muchos factores de \(2\), como \(998244353=119\times 2^{23}+1\).

Código

El código es similar al FFT, con ajustes para aritmética modular.

int n,a[N],ff[N],f[N],g[N];
ll factorial[N], inverse[N], roots[N];
int fastPow(int base, int exp, int mod){
	int result = 1;
	while(exp > 0){
		if(exp & 1) result = (ll)result * base % mod;
		base = (ll)base * base % mod;
		exp >>= 1;
	}
	return result;
}
int modAdd(int x, int mod){
	return x >= mod ? x - mod : x;
}
void nttTransform(int arr[], int len, int dir){
	for(int i=0, j=0; i<len; i++){
		if(i < j) swap(arr[i], arr[j]);
		for(int k=len>>1; (j^=k) < k; k>>=1);
	}
	roots[0] = 1;
	for(int step=1; step<len; step<<=1){
		int step2 = step<<1;
		int base = fastPow(dir==1 ? 3 : 332748118, (mod-1)/step2, mod);
		for(int i=step2-2; i>=0; i-=2){
			roots[i] = roots[i>>1];
			roots[i+1] = roots[i] * base % mod;
		}
		for(int i=0; i<len; i+=step2){
			for(int j=i; j<i+step; j++){
				int val1 = arr[j], val2 = roots[j-i] * arr[j+step] % mod;
				arr[j] = modAdd(val1 + val2, mod);
				arr[j+step] = modAdd(val1 - val2 + mod, mod);
			}
		}
	}
	if(dir == -1){
		int invLen = fastPow(len, mod-2, mod);
		for(int i=0; i<len; i++){
			arr[i] = (ll)arr[i] * invLen % mod;
		}
	}
}

Multiplicación Polinomial con Módulo Arbitrario

Generalmente, el rango de valores es grande, FFT puede tener problemas de precisión, y el módulo no es NTT. Se requieren técnicas especiales.

NTT con Tres Módulos

Se puede usar NTT nueve veces con tres módulos de orden \(10^9\), calculando por separado y combinando con el Teorema Chino del Resto antes de aplicar el módulo arbitrario. Ventaja: conceptualmente simple. Desventaja: alto costo constante.

FFT con Coeficientes Separados

8 FFT

Dividimos cada coeficiente \(a_i\) de \(f,g\) en \(x\times 2^{15}+y\), y luego multiplicamos. Sea \(f=f_1\times 2^{15}+f_2\), similar para \(g\). Entonces \(f*g=f_2*g_2+(f_1*g_2+f_2*g_1)\times 2^{15}+f_1*g_1\times 4^{15}\). Requiere 8 FFT, menos eficiente que 9 NTT.

5 FFT

Aprovechando los números complejos, podemos unir \(f_1,f_2\) en la parte real e imaginaria de un polinomio: \(f'=f_1+f_2I\). Luego, \(f'*g_1\) y \(f'*g_2\) dan los componentes necesarios. Esto requiere 5 FFT, más rápido aún. Existe una variante con 4 FFT, no cubierta aquí.

FWT

Convolución xor/or/and

Sea \(\otimes\in{\operatorname{or},\operatorname{xor},\operatorname{and}}\), calculamos

[C_i=\sum_{j\otimes k=i}A_jB_k ]Análogo a FFT, convertimos a multiplicación punto a punto: \(FWT_C[i]=FWT_A[i]*FWT_B[i]\), luego transformamos inversamente.

Transformaciones específicas:

Convolución or

for(int step=1; step<len; step<<=1){
		int step2 = step<<1;
		for(int i=0; i<len; i+=step2){
			for(int j=i; j<i+step; j++){
				a[j+step] = dir==1 ? modAdd(a[j+step]+a[j], mod) : modAdd(a[j+step]-a[j]+mod, mod);
			}
		}
	}

Convolución and

for(int step=1; step<len; step<<=1){
		int step2 = step<<1;
		for(int i=0; i<len; i+=step2){
			for(int j=i; j<i+step; j++){
				a[j] = dir==1 ? modAdd(a[j]+a[j+step], mod) : modAdd(a[j]-a[j+step]+mod, mod);
			}
		}
	}

Convolución xor

for(int step=1; step<len; step<<=1){
		int step2 = step<<1;
		for(int i=0; i<len; i+=step2){
			for(int j=i; j<i+step; j++){
				int tempX = a[j], tempY = a[j+step];
				a[j] = modAdd(tempX+tempY, mod);
				a[j+step] = modAdd(tempX-tempY+mod, mod);
			}
		}
	}
	if(dir == -1){
		ll invLen = fastPow(len, mod-2, mod);
		for(int i=0; i<len; i++){
			a[i] = invLen * a[i] % mod;
		}
	}

Convolución de Subconjuntos

[C_S=\sum_{T \subseteq S} A_TB_{S-T} ]Se puede reformular como \(C_S=\sum_{i\cap j=\emptyset,i\cup j=S} A_iB_{j}\), similar a la convolución or pero con la condición de intersección vacía. Al expandir, obtenemos \(C_S=\sum_{|i|+|j|=|S|,i\cup j=S} A_iB_{j}\). Entonces, extendemos la convolución or con una dimensión adicional para el tamaño del conjunto, y combinamos.

for(int i=0; i<tamano; i++){
		cin >> a[ __builtin_popcount(i) ][i];
	}
	for(int i=0; i<tamano; i++){
		cin >> b[ __builtin_popcount(i) ][i];
	}
	for(int i=0; i<=n; i++){
		fwtTransform(a[i], tamano, 1);
		fwtTransform(b[i], tamano, 1);
	}
	for(int j=0; j<=n; j++){
		for(int k=0; k<=j; k++){
			for(int i=0; i<tamano; i++){
				c[j][i] = modAdd(c[j][i] + (long long)a[k][i] * b[j-k][i] % mod, mod);
			}
		}
	}
    for(int i=0; i<=n; i++){
		fwtTransform(c[i], tamano, -1);
	}

Operaciones Polinomiales Avanzadas

La mayoría usa NTT para multiplicación.

Inversión Polinomial

Para un polinomio \(F(x)\), encontrar \(G(x)\) tal que \(F(x)*G(x) \equiv 1 \mod x^n\) módulo un primo NTT \(p\).

Enfoque directo

Se puede hacer por recurrencia, complejidad \(\mathcal O(nm)\). La condición necesaria y suficiente es que el término constante de \(F(x)\) tenga inverso módulo \(p\).

Método de duplicación

Usamos duplicación. Sea \(n=2^k\). Asumamos que tenemos \(F(x)*G_m(x)\equiv 1 \mod x^m\). Queremos \(F(x)*G_{2m}(x)\equiv 1 \mod x^{2m}\).

Derivación:

Primero, \(G_{2m}(x)-G_m(x)\equiv 0 \mod x^{m} \). Al cuadrado:

[(G_{2m}(x)-G_m(x))^2\equiv 0 \mod x^{2m} ][G_{2m}(x)^2-2G_{2m}(x)*G_m(x)+G_m(x)^2\equiv 0 \mod x^{2m} ]Usando \(F(x)*G_{2m}(x)\equiv 1 \mod x^{2m} \), multiplicamos por \(F(x)\):

[G_{2m}^2(x)F(x)-2G_{2m}(x)G_m(x)F(x)+G_m(x)^2F(x)\equiv 0 \mod x^{2m} ][G_{2m}(x)-2G_m(x)+G_m(x)^2F(x)\equiv 0 \mod x^{2m} ][G_{2m}(x)\equiv 2G_m(x)-G_m(x)^2F(x) \mod x^{2m} ]El lado derecho es conocido; se calcula con NTT. Complejidad por duplicación: \(\mathcal O(m\log m)\), total \(\mathcal O(n\log n)\).

Optimización: \(G_{2m}(x)\equiv G_m(x)(2-G_m(x)F(x)) \mod x^{2m} \). Notar que \(G_m(x)^2F(x)\) tiene grado hasta \(4m-3\), requiere longitud \(4m\).

Optimización adicional

El método anterior usa 3 NTT de longitud \(4m\). Se puede optimizar. Dividimos \(G_m(x)^2*F(x)\) en dos pasos: primero \(C(x)=F(x)*G_m(x)\mod x^{2m}\), luego \(D(x)=C(x)*G_m(x)\mod x^{2m}\). Observamos que \(C[0]=1\) y \(C[1..m-1]=0\), así que solo necesitamos calcular \(C[m..2m-1]\). Análogo para \(D(x)\). En total, 5 NTT de longitud \(2m\).

División Polinomial

Dados \(F(x)\) de grado \(n\) y \(G(x)\) de grado \(m\), encontrar \(Q(x), R(x)\) con \(Q\) de grado \(n-m\) y \(R\) de grado \([F_R(x)=Q_R(x)*G_R(x)+x^{n-m+1}*R_R(x) ]Modulando por \(x^{n-m+1}\) eliminamos \(R_R(x)\):

[F_R(x)=Q_R(x)*G_R(x)\mod x^{n-m+1} ][Q_R(x)=F_R(x)*G_R^{-1}(x)\mod x^{n-m+1} ]Se calcula la inversa de \(G_R(x)\), se obtiene \(Q_R(x)\), se invierte para \(Q(x)\), y luego \(R(x)=F(x)-G(x)*Q(x)\). Complejidad \(\mathcal O(n\log n)\).

Recurrencia Lineal Homogénea de Coeficientes Constantes

Calcular el \(n\)-ésimo término de una secuencia \(f_i\) que satisface:

[f_n=\sum_{i=1}^k a_if_{k-i} ]### Exponenciación matricial

Directamente con matrices, complejidad \(\mathcal O(k^3\log n)\). Optimización: notar que \(f_j\) se puede expresar como combinación lineal de \(f_1..f_k\). El proceso de expansión se relaciona con reducir potencias de polinomios. Específicamente, \(x^n \mod x^k-\sum_{i=1}^k a_ix^{k-i}\). Se puede usar exponenciación rápida con división polinomial. Si el módulo es NTT, complejidad \(\mathcal O(k\log k\log n)\); si no, con división larga y multiplicación directa, \(\mathcal O(k^2\log n)\).

Raíz Cuadrada Polinomial

Encontrar \(G(x)^2\equiv F(x)\mod x^n\).

Enfoque directo

Similar a la inversión, se puede hacer por recurrencia si \(F[0]\) tiene raíz cuadrada módulo \(p\).

Método de duplicación

Análogo a la inversión. Asumamos \(G_m^2(x)\equiv F(x)\mod x^m\). Queremos \(G_{2m}(x)^2\equiv F(x)\mod x^{2m}\). Como \(G_{2m}(x)-G_m(x)\equiv 0 \mod x^m \), al cuadrado y manipulando:

[G_{2m}(x)=\frac{F(x)+G_m(x)^2}{2G_m(x)}\mod x^{2m} ]Requiere inversión, complejidad \(\mathcal O(n\log n)\).

Evaluación Multi-Punto Polinomial

Para un polinomio \(F(x)\) de grado \(n\), calcular \(F(a_1),F(a_2),\dots,F(a_m)\).

Divide y vencerás

Observar que \(F(a_i)\equiv F(x)\mod x-a_i\). Usamos división recursiva. Para un intervalo \([l,r]\), tenemos \(F'(x)\), calculamos \(F'(x)\mod \prod_{i=l}^r (x-a_i)\) y recursamos. El producto se precalcula con NTT dividido, complejidad \(\mathcal O(n\log^2 n)\). Total \(\mathcal O(n\log^2 n)\).

Logaritmo Polinomial

Repaso de límites, derivadas e integrales...

Dados \(F(x)\), encontrar \(G(x)\equiv \ln F(x)\mod x^n\). Requiere \(F[0]=1\).

Sea \(f(x)=\ln x\), entonces \(G(x)=f(F(x))\). Derivando: \(G'(x)\equiv \frac{F'(x)}{F(x)}\mod x^n\). Calculamos \(G'(x)\) y luego integramos.

Derivada polinomial: \(F'(x)=\sum_{i=0}^{n-1}(i+1)F[i+1]x^i\). Integral: \(\int F(x)=\sum_{i=0}^{n}\frac{F'[i]x^{i+1}}{i+1}\). Entonces \(G(x)=\int \frac{F'(x)}{F(x)}\mod x^n\).

Exponencial Polinomial

Dados \(A(x)\), encontrar \(B(x)\equiv e^{A(x)}\mod x^n\). Requiere \(A[0]=0\).

Toma logaritmo: \(\ln B(x)\equiv A(x)\mod x^n\). Usamos iteración de Newton.

No necesitamos saber la demostración, solo aplicar la fórmula.

Sea \(G(x)=\ln x - A(x)\), donde \(x\) es un polinomio y \(A(x)\) es constante. Aplicamos la fórmula:

[B_{2m}(x)\equiv B_m(x)-\frac{G(B_m(x))}{G'(B_m(x))}\mod x^{2m} ][B_{2m}(x)\equiv B_m(x)-\frac{\ln B_m(x)-A(x)}{1/B_m(x)}\mod x^{2m} ][B_{2m}(x)\equiv B_m(x)-B_m(x)(\ln B_m(x)-A(x))\mod x^{2m} ][B_{2m}(x)\equiv B_m(x)(1-\ln B_m(x)+A(x))\mod x^{2m} ]

Potencia Rápida Polinomial

Dados \(A(x),k\), encontrar \(B(x)\equiv A(x)^k\mod x^n\). \(k\) es un número grande.

Exponenciación binaria directa

Complejidad \(O(n\log n\log k)\), no óptimo si \(k\) es muy grande.

Usando ln y exp

Notar que \(\ln B(x)\equiv k\ln A(x)\mod x^n\), entonces \(B(x)=e^{k\ln A(x)}\). Requiere \(A[0]=1\). Si no es así, ajustamos: sea \(A(x)=px^qC(x)\) con \(C[0]=1\). Entonces \(B(x)\equiv p^kx^{qk}e^{k\ln C(x)}\mod x^n\). Nota: \(p^k\) módulo \(\varphi(p)\), y si \(qk\ge n\), \(B(x)=0\).

Interpolación Rápida

Divide y vencerás

Interpolación de Lagrange rápida. Complejiadd \(\mathcal O(n\log^2 n)\).

La fórmula es:

[F(x)=\sum_{i=1}^{n}y_i\prod_{j\neq i}\frac{x-x_j}{x_i-x_j} ]Reescribimos como:

[F(x)=\sum_{i=1}^{n}\frac{y_i}{\prod_{k\neq i}(x_i-x_k)}\prod_{j\neq i}(x-x_j) ]Sea \(G_i(x)=\prod_{j\neq i}(x-x_j)\), entonces \(\prod_{k\neq i}(x_i-x_k)=G_i(x_i)\). Definimos \(G(x)=\sum_{i=1}^{n}G_i(x)\), que se calcula con NTT dividido. Evaluamos \(G(x_1),\dots,G(x_n)\) con evaluación multi-punto. Finalmente, \(F(x)=\sum_{i=1}^{n}\frac{y_i}{G(x_i)}G_i(x)\), calculable con NTT dividido.

Optimización

El método anterior tiene alta constante. Observando por derivadas: si \(W(x)=\prod_{i=1}^{n}(x-x_i)\), entonces \(W'(x)=\sum_{i=1}^{n}\prod_{j\neq i}(x-x_j)=G(x)\). Esto reduce la constante en aproximadamente un cuarto.

Código

Se incluye código para inversión, división y raíz cuadrada. Se usa vector para operaciones polinomiales, con resize para manejo dinámico.

[Insertar código aquí]