用途在 O ( n l o g n ) O(nlog_n) O(nlogn?)復雜度內解決多項式乘法 比 O ( N 2 ) O(N^2) O(N2)要優
A ( x ) = a 0 + a 1 x + . . . + a n x n A(x)=a_0+a_1x+...+a_nx^n A(x)=a0?+a1?x+...+an?xn
性質一:可以用n+1個點,表示一個n次多項式
證明用高斯消元,范德蒙行列式滿秩唯一解,
點表示法:
如果多項式乘積為: C ( x ) = A ( x ) B ( x ) C(x)=A(x)B(x) C(x)=A(x)B(x)
那么:如果A(x)是n次的,B(x)是m次的,那么我們能用n+m+1個點表示C(x),
x 1 , x 2 , x 3 . . . . . x n + m + 1 x_1,x_2,x_3.....x_{n+m+1} x1?,x2?,x3?.....xn+m+1?
A對應的點為: ( x 1 , A ( x 1 ) ) , ( x 2 , A ( x 2 ) . . . . ( x n + m + 1 , A ( x n + m + 1 ) ) (x_1,A(x_1)),(x_2,A(x_2)....(x_{n+m+1},A(x_{n+m+1})) (x1?,A(x1?)),(x2?,A(x2?)....(xn+m+1?,A(xn+m+1?))
B對應的點為: ( x 1 , B ( x 1 ) ) , ( x 2 , B ( x 2 ) . . . . ( x n + m + 1 , B ( x n + m + 1 ) ) (x_1,B(x_1)),(x_2,B(x_2)....(x_{n+m+1},B(x_{n+m+1})) (x1?,B(x1?)),(x2?,B(x2?)....(xn+m+1?,B(xn+m+1?))
我們知道C(x)=A(x)B(x)
那么C對應的點為: ( x 1 , A ( x 1 ) B ( x 1 ) ) , ( x 2 , A ( x 2 ) B ( x 2 ) . . . . ( x n + m + 1 , A ( x n + m + 1 ) B ( x n + m + 1 ) ) (x_1,A(x_1)B(x_1)),(x_2,A(x_2)B(x_2)....(x_{n+m+1},A(x_{n+m+1})B(x_{n+m+1})) (x1?,A(x1?)B(x1?)),(x2?,A(x2?)B(x2?)....(xn+m+1?,A(xn+m+1?)B(xn+m+1?))
所以我們可以通過O(n+m)表示出來C(x);
所以我們希望從系數表示法轉化成點表示法,在從點表示法轉化成系數表示法,
復數
歐拉公式證明:
證明程序參考繁凡博客吧,
https://www.wolai.com/naS2MSmNf2imHpmYtEpUtF#bedxV3DwQ7wfDDThBQg78S

遞回版:未運行出來:
///FFT遞回版
#include <bits/stdc++.h>
using namespace std;
const double pi = acos(-1);
const int N = 300;
struct Complex
{
double x, y;
Complex(double x = 0, double y = 0): x(x), y(y) {}
} A[N], B[N];
Complex operator * (Complex J, Complex Q)
{
//模長相乘,幅度相加
return Complex(J.x * Q.x - J.y * Q.y, J.x * Q.y + J.y * Q.x);
}
Complex operator - (Complex J, Complex Q)
{
return Complex(J.x - Q.x, J.y - Q.y);
}
Complex operator + (Complex J, Complex Q)
{
return Complex(J.x + Q.x, J.y + Q.y);
}
void FFT(int limit, Complex *a, int type)
{
if(limit == 1) return ;//常數項
Complex a1[limit / 2], a2[limit / 2]; //分成兩端,左右兩端平分
for(int i = 0; i <= limit; i += 2) //偶數項,奇數項
{
a1[i / 2] = a[i];
a2[i / 2] = a[i + 1];
}
FFT(limit / 2, a1, type);
FFT(limit / 2, a2, type);
Complex tmp = Complex(cos(2.0 * pi / limit), type * sin(2.0 * pi / limit)), w = Complex(1, 0);
///tmp表示單位根,w表示冪,
for(int i = 0; i < (limit /2); i++, w = w * tmp)
{
a[i] = a1[i] + w * a2[i]; //左加
a[i + limit / 2] = a1[i] - w * a2[i]; //右減
}
}
int main()
{
int n, m;
scanf("%d%d", &n, &m);
for(int i = 0; i <= n; i++)scanf("%lf", &A[i].x);
for(int i = 0; i <= m; i++)scanf("%lf", &B[i].x);
int limit = 1;
while(limit <= n + m) limit *= 2; //2的整數次冪
FFT(limit, A, 1);
FFT(limit, B, 1);
//后面的1表示要進行的變換是什么型別
// 1表示從系數變為點值
//-1表示從點值變為系數
//這就與推導程序中的 w指數部分正負有關,
for(int i = 0; i <= limit; i++)
{
A[i] = A[i] * B[i];
}
FFT(limit, A, -1);
for(int i = 0;i <= n + m; i++)
{
printf("%d ", (int)(A[i].x / limit + 0.5)); ///對其四舍五入+0.5
}
return 0;
}
迭代版
蝴蝶變換
///迭代版
#include <bits/stdc++.h>
using namespace std;
const int N = 5e6 + 7;
const double PI = acos(-1);
int n, m;
int limit = 1;
int res, ans[N];
int l;
int r[N];
struct Complex
{
double x, y;
Complex (double x = 0, double y = 0) : x(x), y(y) {}
} a[N], b[N];
Complex operator * (Complex J, Complex Q)
{
//模長相乘,幅度相加
return Complex(J.x * Q.x - J.y * Q.y, J.x * Q.y + J.y * Q.x);
}
Complex operator + (Complex J, Complex Q)
{
return Complex(J.x + Q.x, J.y + Q.y);
}
Complex operator - (Complex J, Complex Q)
{
return Complex(J.x - Q.x, J.y - Q.y);
}
void FFT(Complex *A, int type)
{
for(int i = 0; i < limit; i++)
{
if(i < r[i]) swap(A[i], A[r[i]]); ///保證不在換過去
}
///從底層開始合并:
for(int mid = 1; mid < limit; mid = mid * 2)
{
///待合并區間長度的一半,最開始是長度為1的合并,類似倍增的思想
Complex wn(cos(PI / mid), type * sin(PI / mid)); //單位根
for(int len = mid *2, pos = 0; pos < limit; pos += len)
{
///len是區間的長度,pos是當前的位置
Complex w(1, 0);
for(int k = 0; k < mid; k++, w = w * wn)
{
///只掃左半部分,蝴蝶變換得到有半部分,
Complex x = A[pos + k]; //左半部分
Complex y = w * A[pos + mid + k]; //右半部分
A[pos + k] = x + y;//左加
A[pos + mid + k] = x - y;///右減
}
}
}
if(type == 1) return ;
for(int i = 0; i <= limit; i++)
{
A[i].x = A[i].x / limit;
A[i].y = A[i].y / limit;///這個版本沒用,加不加沒影響,優化的版本必須要加
///除以我們推出的N,
}
}
int main()
{
cin >> n >> m;
for(int i = 0; i <= n; i++) scanf("%lf", &a[i].x);
for(int j = 0; j <= m; j++) scanf("%lf", &b[j].x);
while(limit <= n + m) limit <<= 1, l++;///求出位數l,和2的整數次冪
for(int i = 0; i < limit; i++)
{
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));///蝴蝶變換,
}
FFT(a, 1);
FFT(b, 1);
for(int i = 0; i <= limit; ++i)
{
a[i] = a[i] * b[i];
}
FFT(a, -1);
for(int i = 0; i <= m + n; i++)
{
printf("%d ", int(a[i].x + 0.5));
}
return 0;
}
優化:三步變兩步:
設a和b是實多項式, F = a + b i F=a+bi F=a+bi則 F 2 = a 2 ? b 2 + 2 a b i F^2=a^2-b^2+2abi F2=a2?b2+2abi ,注意我們要求的是F虛部的一半,這樣我們通過兩次FFT就可以求出結果了,
所以我們可以把 b ( x ) b(x) b(x)放到a(x)的虛部,求出 a ( x ) 2 a(x)^2 a(x)2然后取出a(x)虛部除以2就是答案,
///優化版
#include <bits/stdc++.h>
using namespace std;
const int N = 5e6 + 7;
const double PI = acos(-1);
int n, m;
int limit = 1;
int res, ans[N];
int l;
int r[N];
struct Complex
{
double x, y;
Complex (double x = 0, double y = 0) : x(x), y(y) {}
} a[N], b[N];
Complex operator * (Complex J, Complex Q)
{
//模長相乘,幅度相加
return Complex(J.x * Q.x - J.y * Q.y, J.x * Q.y + J.y * Q.x);
}
Complex operator + (Complex J, Complex Q)
{
return Complex(J.x + Q.x, J.y + Q.y);
}
Complex operator - (Complex J, Complex Q)
{
return Complex(J.x - Q.x, J.y - Q.y);
}
void FFT(Complex *A, int type)
{
for(int i = 0; i < limit; i++)
{
if(i < r[i]) swap(A[i], A[r[i]]); ///保證不在換過去
}
///從底層開始合并:
for(int mid = 1; mid < limit; mid = mid * 2)
{
///待合并區間長度的一半,最開始是長度為1的合并,類似倍增的思想
Complex wn(cos(PI / mid), type * sin(PI / mid)); //單位根
for(int len = mid *2, pos = 0; pos < limit; pos += len)
{
///len是區間的長度,pos是當前的位置
Complex w(1, 0);
for(int k = 0; k < mid; k++, w = w * wn)
{
///只掃左半部分,蝴蝶變換得到有半部分,
Complex x = A[pos + k]; //左半部分
Complex y = w * A[pos + mid + k]; //有半部分
A[pos + k] = x + y;
A[pos + mid + k] = x - y;
}
}
}
if(type == 1) return ;
for(int i = 0; i <= limit; i++)
{
A[i].x = A[i].x / limit;
A[i].y = A[i].y / limit;
///除以我們推出的N,
}
}
int main()
{
cin >> n >> m;
for(int i = 0; i <= n; i++) scanf("%lf", &a[i].x);
for(int j = 0; j <= m; j++) scanf("%lf", &a[j].y);
while(limit <= n + m) limit <<= 1, l++;
for(int i = 0; i < limit; i++)
{
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
}
FFT(a, 1);
///FFT(b, 1);
for(int i = 0; i <= limit; ++i)
{
a[i] = a[i] * a[i];
}
FFT(a, -1);
for(int i = 0; i <= m + n; i++)
{
printf("%d ", int(a[i].y/2 + 0.5));
}
return 0;
}
例題:P1919 【模板】FFT快速傅里葉變換
///優化版
#include <bits/stdc++.h>
using namespace std;
const int N = 5e6 + 7;
const double PI = acos(-1);
int n, m;
int limit = 1;
int res, ans[N];
int l;
int r[N];
struct Complex
{
double x, y;
Complex (double x = 0, double y = 0) : x(x), y(y) {}
} a[N], b[N];
Complex operator * (Complex J, Complex Q)
{
//模長相乘,幅度相加
return Complex(J.x * Q.x - J.y * Q.y, J.x * Q.y + J.y * Q.x);
}
Complex operator + (Complex J, Complex Q)
{
return Complex(J.x + Q.x, J.y + Q.y);
}
Complex operator - (Complex J, Complex Q)
{
return Complex(J.x - Q.x, J.y - Q.y);
}
void FFT(Complex *A, int type)
{
for(int i = 0; i < limit; i++)
{
if(i < r[i]) swap(A[i], A[r[i]]); ///保證不在換過去
}
///從底層開始合并:
for(int mid = 1; mid < limit; mid = mid * 2)
{
///待合并區間長度的一半,最開始是長度為1的合并,類似倍增的思想
Complex wn(cos(PI / mid), type * sin(PI / mid)); //單位根
for(int len = mid * 2, pos = 0; pos < limit; pos += len)
{
///len是區間的長度,pos是當前的位置
Complex w(1, 0);
for(int k = 0; k < mid; k++, w = w * wn)
{
///只掃左半部分,蝴蝶變換得到有半部分,
Complex x = A[pos + k]; //左半部分
Complex y = w * A[pos + mid + k]; //有半部分
A[pos + k] = x + y;
A[pos + mid + k] = x - y;
}
}
}
if(type == 1) return ;
for(int i = 0; i <= limit; i++)
{
ans[i] += (int)(A[i].y / limit/2 + 0.5);
if(ans[i] >= 10) ///考慮進位問題,
{
ans[i + 1] = ans[i] / 10;
ans[i] %= 10;
limit += (i == limit); ///進位
}
///除以我們推出的N,
}
}
string str, s;
int main()
{
cin >> str >> s;
n = str.size();
m = s.size();
int cnt, tot;
cnt = tot = 0;
for(int i = n - 1; i >= 0; i--) a[cnt++].x = str[i] - '0';
for(int i = m - 1; i >= 0; i--) a[tot++].y = s[i] - '0';
while(limit < n + m) limit <<= 1, l++;
for(int i = 0; i <= limit; i++)
{
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
}
FFT(a, 1);
///FFT(b, 1);
for(int i = 0; i <= limit; ++i)
{
a[i] = a[i] * a[i];
}
FFT(a, -1);
while(!ans[limit] && limit >= 1) limit--; ///去除前導0
limit++;
while(--limit >= 0)
{
cout << ans[limit];
}
return 0;
}
正常版,
//正常的
#include <bits/stdc++.h>
using namespace std;
const int N = 5e6 + 7;
const double PI = acos(-1);
int n, m;
int limit = 1;
int res, ans[N];
int l;
int r[N];
struct Complex
{
double x, y;
Complex (double x = 0, double y = 0) : x(x), y(y) {}
} a[N], b[N];
Complex operator * (Complex J, Complex Q)
{
//模長相乘,幅度相加
return Complex(J.x * Q.x - J.y * Q.y, J.x * Q.y + J.y * Q.x);
}
Complex operator + (Complex J, Complex Q)
{
return Complex(J.x + Q.x, J.y + Q.y);
}
Complex operator - (Complex J, Complex Q)
{
return Complex(J.x - Q.x, J.y - Q.y);
}
void FFT(Complex *A, int type)
{
for(int i = 0; i < limit; i++)
{
if(i < r[i]) swap(A[i], A[r[i]]); ///保證不在換過去
}
///從底層開始合并:
for(int mid = 1; mid < limit; mid = mid * 2)
{
///待合并區間長度的一半,最開始是長度為1的合并,類似倍增的思想
Complex wn(cos(PI / mid), type * sin(PI / mid)); //單位根
for(int len = mid * 2, pos = 0; pos < limit; pos += len)
{
///len是區間的長度,pos是當前的位置
Complex w(1, 0);
for(int k = 0; k < mid; k++, w = w * wn)
{
///只掃左半部分,蝴蝶變換得到有半部分,
Complex x = A[pos + k]; //左半部分
Complex y = w * A[pos + mid + k]; //有半部分
A[pos + k] = x + y;
A[pos + mid + k] = x - y;
}
}
}
if(type == 1) return ;
for(int i = 0; i <= limit; i++)
{
ans[i] += (int)(A[i].x / limit + 0.5);
if(ans[i] >= 10) ///考慮進位問題,
{
ans[i + 1] = ans[i] / 10;
ans[i] %= 10;
limit += (i == limit); ///進位
}
///除以我們推出的N,
}
}
string str, s;
int main()
{
cin >> str >> s;
n = str.size();
m = s.size();
int cnt, tot;
cnt = tot = 0;
for(int i = n - 1; i >= 0; i--) a[cnt++].x = str[i] - '0';
for(int i = m - 1; i >= 0; i--) b[tot++].x = s[i] - '0';
while(limit < n + m) limit <<= 1, l++;
for(int i = 0; i <= limit; i++)
{
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
}
FFT(a, 1);
FFT(b, 1);
for(int i = 0; i <= limit; ++i)
{
a[i] = a[i] * b[i];
}
FFT(a, -1);
while(!ans[limit] && limit >= 1) limit--; ///去除前導0
limit++;
while(--limit>=0)
{
cout << ans[limit];
}
return 0;
}
P3338 [ZJOI2014]力
先放上,防止自己咕咕咕
題意:
n個數
q 1 , q 2 , . . . q n q_1,q_2,...q_n q1?,q2?,...qn?定義
F j = ∑ i = 1 j ? 1 q i × q j ( i ? j ) 2 ? ∑ i = j + 1 n q i × q j ( i ? j ) 2 F_j = \sum_{i=1}^{j?1}\frac{q_i×q_j}{(i?j)^2} ? \sum_{i=j+1}^n\frac{q_i×q_j}{(i?j)^2} Fj?=∑i=1j?1?(i?j)2qi?×qj???∑i=j+1n?(i?j)2qi?×qj??
E i = F i q i E_i~=~\frac{F_i}{q_i} Ei? = qi?Fi??
讓求 E i E_i Ei?




快速數論變換 (NTT)
#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int N = 5e6 + 7;
const double PI = acos(-1);
const int p = 998244353, G = 3, Gi = 332748118;
int n, m;
int limit = 1;
int res, ans[N];
int l;
int r[N];
ll a[N], b[N];
ll qpow(ll a, ll b)
{
ll ans = 1;
while(b)
{
if(b & 1) ans = ans * a % p;
a = a * a % p;
b >>= 1;
}
return ans;
}
void NTT(ll *A, int type)
{
for(int i = 0; i < limit; i++)
{
if(i < r[i]) swap(A[i], A[r[i]]); ///保證不在換過去
}
///從底層開始合并:
for(int mid = 1; mid < limit; mid = mid * 2)
{
///待合并區間長度的一半,最開始是長度為1的合并,類似倍增的思想
ll wn = qpow(G, (p - 1) / (mid * 2));
if(type==-1) wn = qpow(wn, p - 2);
for(int len = mid * 2, pos = 0; pos < limit; pos += len)
{
///len是區間的長度,pos是當前的位置
ll w = 1;
for(int k = 0; k < mid; k++, w = w * wn % p)
{
///只掃左半部分,蝴蝶變換得到有半部分,
int x = A[pos + k]; //左半部分
int y = w * A[pos + mid + k] % p; //有半部分
A[pos + k] = (x + y) % p;
A[pos + mid + k] = (x - y + p) % p;
}
}
}
if(type == 1) return ;
ll inlimit = qpow(limit, p - 2);
for(int i = 0; i < limit; i++)
{
A[i] = (A[i] * inlimit) % p;
///除以我們推出的N,
}
}
int main()
{
cin >> n >> m;
for(int i = 0; i <= n; i++) cin >> a[i], a[i] = (a[i] + p) % p;
for(int i = 0; i <= m; i++) cin >> b[i], b[i] = (b[i] + p) % p;
while(limit <= n + m) limit <<= 1, l++;
for(int i = 0; i <= limit; i++)
{
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
}
NTT(a, 1);
NTT(b, 1);
for(int i = 0; i <= limit; ++i)
{
a[i] = a[i] * b[i] % p;
}
NTT(a, -1);
for(int i = 0; i <= n + m; i++) cout << a[i] << " ";
return 0;
}
轉載請註明出處,本文鏈接:https://www.uj5u.com/qita/298371.html
標籤:其他
上一篇:LeetCode二叉樹的層序遍歷
