1. 程式人生 > 其它 >FFT,FWT演算法比對——HDU6966. I love sequences(FWT)

FFT,FWT演算法比對——HDU6966. I love sequences(FWT)

卷積定義:現有\(A,B\)兩個序列,\(A*B\)為其卷積

\[(A*B)_k=\sum\limits_{i+j=k}A_iB_j \]

對其進行拓展,可以將\(i+j=k\)中的加號換成其它的二元運算子,例如與,或,異或等等。

FFT——快速傅立葉變換

快速傅立葉變換就是針對\(i+j=k\)型別的卷積的。

\(T\)為傅立葉變換,我們要求這個變換滿足

\[T(A)_kT(B)_k=T(A*B)_k \]

這樣計算\(A*B\)就可以通過計算\(T(A),T(B)\)得到\(T(A*B)\),再對其進行\(T^{-1}\)變換得到。

\(T\)是一種線性變換,所以我們可以將\(T(A)_k\)

表示成如下的式子

\[T(A)_k=\sum\limits_{i=0}^{n-1}A_if(k,i) \]

要使得\(T(A)_kT(B)_k=T(A*B)_k\)成立,就要

\[\begin{aligned} \left(\sum\limits_{i=0}^{n-1}A_if(k,i)\right)\left(\sum\limits_{j=0}^{n-1}B_jf(k,j)\right)&=\sum\limits_{m=0}^{n-1}\sum\limits_{i+j=m}A_iB_jf(k,m)\\ \sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}A_iB_jf(k,i)f(k,j)&=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}A_iB_jf(k,i+j)\\ \end{aligned} \]

即需要

\[f(k,i)f(k,j)=f(k,i+j) \]

顯然我們可以令\(f(k,i)=a^i\),理論上\(a\)的取值可以任意。

FFT中,選擇了\(n\)次單位根的冪次\(\omega_n^m(0\leq m<n)\)作為\(a\)

這樣保證了有\(n\)組相互之間線性無關的線性變換,這保證了逆變換還能變得回來。

FFT採用分治演算法進行,根據奇偶,將原序列分成兩部分,\(A_{(\dots 0)_2}\)以及\(A_{(\dots 1)_2}\)

\[\begin{aligned} T(A)_k&=\sum\limits_{i=0}^{n-1}A_i\omega_n^{ki}\\ &=\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i}\omega_n^{2ki}+\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i+1}\omega_n^{2ki+k}\\ &=\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i}\omega_{\frac{n}{2}}^{ki}+\omega_n^k\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i+1}\omega_\frac{n}{2}^{ki}\\ &=T(A_{(\dots 0)_2})_k+\omega_n^kT(A_{(\dots 1)_2})_k&(k<\frac{n}{2})\\ T(A)_k&=\sum\limits_{i=0}^{n-1}A_i\omega_n^{ki}\\ &=\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i}\omega_n^{2ki}+\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i+1}\omega_n^{2ki+k}\\ &=\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i}\omega_{\frac{n}{2}}^{ki}+\omega_n^k\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i+1}\omega_\frac{n}{2}^{ki}\\ &=\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i}\omega_{\frac{n}{2}}^{(k-\frac{n}{2})i}+\omega_n^k\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i+1}\omega_\frac{n}{2}^{(k-\frac{n}{2})i}\\ &=T(A_{(\dots 0)_2})_{k-\frac{n}{2}}+\omega_n^kT(A_{(\dots 1)_2})_{k-\frac{n}{2}}\\ &=T(A_{(\dots 0)_2})_{k-\frac{n}{2}}-\omega_n^{k-\frac{n}{2}}T(A_{(\dots 1)_2})_{k-\frac{n}{2}}&(k\geq\frac{n}{2})\\ \end{aligned} \]

這樣就可以進行FFT

了,逆變換求逆矩陣。

最終形成了一個線性變換,

不難得出,最終的變換矩陣為

\[A=\begin{bmatrix} (\omega_n^0)^0&(\omega_n^0)^1&\cdots&(\omega_n^0)^{n-1}\\ (\omega_n^1)^0&(\omega_n^1)^1&\cdots&(\omega_n^1)^{n-1}\\ \vdots&\vdots&\ddots&\vdots\\ (\omega_n^{n-1})^0&(\omega_n^{n-1})^1&\cdots&(\omega_n^{n-1})^{n-1} \end{bmatrix}\\ \]

其逆矩陣

\[A^{-1}=\frac{1}{n}\begin{bmatrix} (\omega_n^0)^0&(\omega_n^0)^1&\cdots&(\omega_n^0)^{n-1}\\ (\omega_n^{-1})^0&(\omega_n^{-1})^1&\cdots&(\omega_n^{-1})^{n-1}\\ \vdots&\vdots&\ddots&\vdots\\ (\omega_n^{1-n})^0&(\omega_n^{1-n})^1&\cdots&(\omega_n^{1-n})^{n-1} \end{bmatrix}\\ \]

FWT——快速沃爾什變換

快速沃爾什變換是針對\(i\oplus j=k\),這裡的\(\oplus\)不一定指異或,可以是任意的二元按位運算子(不一定是二進位制,但下面以二進位制為例)

\(T\)為沃爾什變換,我們要求這個變換滿足

\[T(A)_kT(B)_k=T(A*B)_k \]

這樣計算\(A*B\)就可以通過計算\(T(A),T(B)\)得到\(T(A*B)\),再對其進行\(T^{-1}\)變換得到。

\(T\)是一種線性變換,所以我們可以將\(T(A)_k\)表示成如下的式子

\[T(A)_k=\sum\limits_{i=0}^{n-1}A_if(k,i) \]

要使得\(T(A)_kT(B)_k=T(A*B)_k\)成立,就要

\[\begin{aligned} \left(\sum\limits_{i=0}^{n-1}A_if(k,i)\right)\left(\sum\limits_{j=0}^{n-1}B_jf(k,j)\right)&=\sum\limits_{m=0}^{n-1}\sum\limits_{i\oplus j=m}A_iB_jf(k,m)\\ \sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}A_iB_jf(k,i)f(k,j)&=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}A_iB_jf(k,i\oplus j)\\ \end{aligned} \]

即需要

\[f(k,i)f(k,j)=f(k,i\oplus j) \]

這裡以二進位制異或為例,認為\(\oplus\)為二進位制異或。

由於是按位運算,我們不妨逐位考慮,令\(f(k,i)=\prod\limits_x f(k,i_x)\)

那麼\(f(k,i)f(k,j)=f(k,i\oplus j)\)成立當且僅當下面\(3\)個式子成立

\[f(k,0)f(k,0)=f(k,0)\\ f(k,1)f(k,1)=f(k,0)\\ f(k,1)f(k,0)=f(k,1) \]

可以解得這\(3\)組解(並不是任取一組就好的,具體原因下面說)

\[f(k,0)=0,f(k,1)=0\\ f(k,0)=1,f(k,1)=-1\\ f(k,0)=1,f(k,1)=1 \]

同樣地,可以採用分治的辦法,根據最高位可以是\(0\)\(1\)可以將序列分為\(A_{(0\dots)_2},A_{(1\dots)_2}\)

\[\text{$i_0$為$i$的最高位,$i_1$為$i$去掉最高位}\\ \begin{aligned} T(A)_k&=\sum\limits_{i=0}^{n-1}A_if(k,i)\\ &=\sum\limits_{i=0}^{\frac{n}{2}-1}A_if(k,i)+\sum\limits_{i=\frac{n}{2}}^{n-1}A_if(k,i)\\ &=f(k_0,0)\sum\limits_{i=0}^{\frac{n}{2}-1}A_if(k_1,i_1)+f(k_0,1)\sum\limits_{i=0}^{\frac{n}{2}-1}A_if(k_1,i_1)\\ &=f(k_0,0)T(A_{(0\dots)_2})_{k_1}+f(k_0,1)T(A_{(1\dots)_2})_{k_1} \end{aligned} \]

注意,式子形式雖然一樣,但是對於\(k<\frac{n}{2}\)\(k\geq\frac{n}{2}\)兩種情況,如果\(f(k_0,0),f(k_0,1)\)取的值分別相同,那麼會造成前一半和後一半是一樣的,顯然變換無法保持可逆。

所以我們要在上面解出的若干個解中,選出\(2\)個線性無關的解,這裡別無選擇,只能選擇

\[f(k,0)=1,f(k,1)=-1\\ f(k,0)=1,f(k,1)=1 \]

讓其中之一(隨便哪一個)作為\(k<\frac{n}{2}\)時的\(f(k_0,0),f(k_0,1)\),另外一組作為\(k\geq\frac{n}{2}\)時的\(f(k_0,0),f(k_0,1)\)

最終線性變換不太好寫,可以寫每一步的矩陣

\[A=\begin{bmatrix} 1&-1\\ 1&1 \end{bmatrix} \]

逆矩陣

\[A^{-1}=\begin{bmatrix} \frac{1}{2}&\frac{1}{2}\\ -\frac{1}{2}&\frac{1}{2} \end{bmatrix} \]

HDU6966. I love sequences

題面:

問題描述:

給定三個序列\(a=[a_1,a_2,\dots,a_n],b=[b_1,b_2,\dots,b_n],c=[c_1,c_2,\dots,c_n]\)每一個含有\(n\)個非負整數。你需要計算

\[\sum\limits_{p=1}^{n}\sum\limits_{k=1}^{+\infty}d_{p,k}c_p^k \]

其中序列\(d=[d_{p,1},d_{p,2},\dots,d_{p,n}]\)按以下的規則計算:

對於每一個整數\(p(1\leq p\leq n)\),每一個整數\(k(1\leq k\leq n)\),有

\[d_{p,k}=\sum\limits_{k=i\oplus j}a_ib_j\\ \text{$i,j$的範圍為$1\leq i\leq \frac{n}{p},1\leq j\leq \frac{n}{p}$} \]

其中\(\oplus\)代表一個三進位制二元按位運算子,滿足\(k_x=\gcd(i_x,j_x)\)

比如\(12\oplus 7=(110)_3\oplus (021)_3=(111)_3=13\)

輸入:

這個問題只有一個樣例。

第一行包含一個整數\(n(n\leq2\times10^5)\),表示序列\(a,b,c\)的長度

第二行包含\(n\)個整數,表示\(a_1,a_2,\dots,a_n\)

第三行包含\(n\)個整數,表示\(b_1,b_2,\dots,b_n\)

第四行包含\(n\)個整數,表示\(c_1,c_2,\dots,c_n\)

輸出:

輸出一個整數,為答案模\(10^9+7\)

樣例輸入:

4
2 3 4 5
1 2 3 4
9 8 7 6

樣例輸出:

1643486

分析:

位運算卷積考慮FWT,考慮\(f(k,i)f(k,j)=f(k,i\oplus j)\)成立的充要條件

\[f(k,0)f(k,0)=f(k,0)\\ f(k,1)f(k,1)=f(k,1)\\ f(k,2)f(k,2)=f(k,2)\\ f(k,0)f(k,1)=f(k,1)\\ f(k,0)f(k,2)=f(k,2)\\ f(k,1)f(k,2)=f(k,1)\\ \]

解得

\[f(k,0)=0,f(k,1)=0,f(k,2)=0\\ f(k,0)=1,f(k,1)=0,f(k,2)=0\\ f(k,0)=1,f(k,1)=0,f(k,2)=1\\ f(k,0)=1,f(k,1)=1,f(k,2)=1\\ \]

同樣地,使用分治

\[\text{$i_0$為$i$的最高位,$i_1$為$i$去掉最高位}\\ \begin{aligned} T(A)_k&=\sum\limits_{i=0}^{n-1}A_if(k,i)\\ &=\sum\limits_{i=0}^{\frac{n}{3}-1}A_if(k,i)+\sum\limits_{i=\frac{n}{3}}^{\frac{2n}{3}-1}A_if(k,i)+\sum\limits_{i=\frac{2n}{3}}^{n-1}A_if(k,i)\\ &=f(k_0,0)\sum\limits_{i=0}^{\frac{n}{3}-1}A_if(k_1,i_1)+f(k_0,1)\sum\limits_{i=\frac{n}{3}}^{\frac{2n}{3}-1}A_if(k_1,i_1)+f(k_0,2)\sum\limits_{i=\frac{2n}{3}}^{n-1}A_if(k_1,i_1)\\ &=f(k_0,0)T(A_{(0\dots)_3})_{k_1}+f(k_0,1)T(A_{(1\dots)_3})_{k_1}+f(k_0,2)T(A_{(2\dots)_3})_{k_1} \end{aligned} \]

\(f(k_0,x)\)需要根據\(k<\frac{n}{3},\frac{n}{3}\leq k<\frac{2n}{3},\frac{2n}{3}\leq k<n\)三種情況,選擇三個線性無關的取值。

所以只能取

\[f(k,0)=1,f(k,1)=0,f(k,2)=0\\ f(k,0)=1,f(k,1)=0,f(k,2)=1\\ f(k,0)=1,f(k,1)=1,f(k,2)=1\\ \]

形成矩陣

\[A=\begin{bmatrix} 1&0&0\\ 1&0&1\\ 1&1&1 \end{bmatrix} \]

其逆矩陣

\[A^{-1}=\begin{bmatrix} 1&0&0\\ 0&-1&1\\ -1&1&0 \end{bmatrix} \]

程式碼:

#include <cstdio>
typedef long long Lint;
const int maxn = 540000;
const int mod = 1e9 + 7;
int a[maxn], b[maxn], c[maxn];
int aa[maxn], bb[maxn], cc[maxn];
void FWT(int* a, int n, int flag) {
    for (int len = 1; len < n; len *= 3) {
        for (int i = 0; i < n; i += len * 3) {
            for (int j = 0; j < len; j++) {
                int x = a[i + j], y = a[i + j + len], z = a[i + j + 2 * len];
                if (flag == 1) {
                    a[i + j + len] = (x + z) % mod;
                    a[i + j + 2 * len] = ((x + y) % mod + z) % mod;
                } else {
                    a[i + j + len] = (mod - y + z) % mod;
                    a[i + j + 2 * len] = (mod - x + y) % mod;
                }
            }
        }
    }
}
int cal(int n, int c) {
    int tot = 1;
    while (tot <= n)
        tot *= 3;
    for (int i = 1; i <= n; i++) {
        aa[i] = a[i];
        bb[i] = b[i];
    }
    for (int i = n + 1; i < tot; i++) {
        aa[i] = bb[i] = 0;
    }
    FWT(aa, tot, 1), FWT(bb, tot, 1);
    for (int i = 0; i < tot; i++) {
        cc[i] = (Lint)aa[i] * bb[i] % mod;
    }
    FWT(cc, tot, -1);
    int res = 0;
    for (int k = 1, ck = c; k < tot; k++, ck = (Lint)ck * c % mod) {
        res = (res + (Lint)cc[k] * ck % mod) % mod;
    }
    return res;
}
int main() {
    int n;
    scanf("%d", &n);
    for (int i = 1; i <= n; i++) {
        scanf("%d", a + i);
    }
    for (int i = 1; i <= n; i++) {
        scanf("%d", b + i);
    }
    for (int i = 1; i <= n; i++) {
        scanf("%d", c + i);
    }
    int ans = 0;
    for (int p = 1; p <= n; p++) {
        ans = (ans + cal(n / p, c[p])) % mod;
    }
    printf("%d\n", ans);
    return 0;
}