2022.03.01 FFT与NNT
https://www.luogu.com.cn/blog/wangrx/solution-p1919
https://www.luogu.com.cn/blog/attack/solution-p38032
P3803 【模板】多项式乘法（FFT）（NTT模板）
https://www.luogu.com.cn/problem/P3803
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=3e6+10;
const int mod=998244353;
const int G=3;
const int Gi=332748118;
int n,m,limit=1,len,rev[N],a[N],b[N];
inline int read(){
int s=0,w=1;
char ch=getchar();
while(ch<'0'||ch>'9'){
if(ch=='-')w=-1;
ch=getchar();
}
while(ch<='9'&&ch>='0'){
s=s*10+ch-'0';
ch=getchar();
}
return s*w;
}
inline int powi(int x,int y){
int fin=1;
while(y){
if(y&1)fin=fin*x%mod;
x=x*x%mod;
y>>=1;
}
return fin;
}
inline void NTT(int *A,int flag){
for(int i=0;i<limit;i++)if(i<rev[i])swap(A[i],A[rev[i]]);
for(int mid=1;mid<limit;mid<<=1){
int Wn=powi(flag==1?G:Gi,(mod-1)/(mid<<1));
for(int j=0;j<limit;j+=(mid<<1)){
int W=1;
for(int k=0;k<mid;k++,W=W*Wn%mod){
int x=A[j+k],y=W*A[j+k+mid]%mod;
A[j+k]=(x+y)%mod;
A[j+k+mid]=(x+mod-y)%mod;
}
}
}
}
signed main(){
n=read();m=read();
for(int i=0;i<=n;i++)a[i]=read(),a[i]=(a[i]+mod)%mod;
for(int i=0;i<=m;i++)b[i]=read(),b[i]=(b[i]+mod)%mod;
while(limit<=n+m)limit<<=1,++len;
for(int i=0;i<limit;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
NTT(a,1);NTT(b,1);
for(int i=0;i<limit;i++)a[i]=a[i]*b[i]%mod;
NTT(a,-1);
int inv=powi(limit,mod-2);
for(int i=0;i<=n+m;i++)cout<<a[i]*inv%mod<<" ";
return 0;
}
P1919 【模板】A*B Problem 升级版（FFT 快速傅里叶变换（注意：可能结果的第一个系数也要进一……）（NTT模板1.1）
https://www.luogu.com.cn/problem/P1919
这里对于NTT模板进行了优化：
- 不用每次计算
powi(flag==1?3:inv3,(mod-1)/(mid<<1))了，直接搞了两个数组，一个叫powi3就是存3的次方的，一个叫powinv3就是存3的次方的逆元的 - 没了/斜眼笑
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=4e6+10;
const int mod=998244353;
int n,a[N],b[N],c[N],rev[N],lena,lenb;
int pow3[N],powinv3[N],inv3,inv,limit=1,len;
char s[N];
inline int powi(int x,int y){
int fin=1;
while(y){
if(y&1)fin=fin*x%mod;
x=x*x%mod;y>>=1;
}
return fin;
}
inline void NTT(int *A,int flag){
for(int i=0;i<limit;i++)if(i<rev[i])swap(A[i],A[rev[i]]);
for(int mid=1;mid<limit;mid<<=1){
int Wn=(flag==1?pow3[mid<<1]:powinv3[mid<<1]);
//int Wn=powi(flag==1?3:inv3,(mod-1)/(mid<<1));
for(int j=0;j<limit;j+=(mid<<1)){
int W=1;
for(int k=0;k<mid;k++,W=(W*Wn)%mod){
int x=A[k+j],y=A[k+j+mid]*W%mod;
A[j+k]=(x+y)%mod;
A[j+k+mid]=(x-y+mod)%mod;
}
}
}
}
signed main(){
inv3=powi(3,mod-2);
//cout<<inv3<<endl;
scanf("%s",s);//printf("%s\n",s);
lena=strlen(s);
for(int i=0;i<lena;i++)a[i]=s[i]-'0';//,cout<<a[i]<<" ";cout<<endl;
scanf("%s",s);//printf("%s\n",s);
lenb=strlen(s);
for(int i=0;i<lenb;i++)b[i]=s[i]-'0';//,cout<<b[i]<<" ";cout<<endl;
--lena;--lenb;
while(limit<=lena+lenb)limit<<=1,++len;
for(int i=0;i<limit;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
for(int i=1;i<=limit;i<<=1)
pow3[i]=powi(3,(mod-1)/i),powinv3[i]=powi(inv3,(mod-1)/i);
NTT(a,1);NTT(b,1);
for(int i=0;i<limit;i++)a[i]=a[i]*b[i]%mod;
NTT(a,-1);
inv=powi(limit,mod-2);
for(int i=0;i<=lena+lenb;i++)a[i]=a[i]*inv%mod;//,cout<<a[i]<<" ";cout<<endl;
for(int i=0;i<=lena+lenb;i++)c[i]=a[lena+lenb-i];//,cout<<c[i]<<" ";cout<<endl;
for(int i=0;i<=lena+lenb+2;i++)if(c[i]>=10)c[i+1]+=c[i]/10,c[i]%=10;
n=lena+lenb+4;
while(c[n]==0)--n;
for(int i=n;i>=0;i--)cout<<c[i];
return 0;
}
P2553 [AHOI2001]多项式乘法
https://www.luogu.com.cn/problem/P2553
这是道美丽的字符串，但是eleveni挂在NTT上
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=666;
const int mod=998244353;
int limit=1,a[N],b[N],rev[N],len,inv,inv3,powi3[N],powinv3[N];
struct node{
int x,y;
bool operator <(const node &b)const{
return y<b.y;
}
}s[N];
string t;
inline int powi(int x,int y){
int fin=1;
while(y){
if(y&1)fin=fin*x%mod;
x=x*x%mod;y>>=1;
}
return fin;
}
inline void NTT(int *A,int flag){
for(int i=0;i<limit;i++)if(i<rev[i])swap(A[i],A[rev[i]]);
//cout<<"Case mid "<<endl;
for(int mid=1;mid<limit;mid<<=1){
//int Wn=flag==1?powi3[mid<<1]:powinv3[mid<<1];
int Wn=powi(flag==1?3:inv3,(mod-1)/(mid<<1));
for(int j=0;j<limit;j+=(mid<<1)){
int W=1;
for(int k=0;k<mid;k++,W=W*Wn%mod){
int x=A[j+k],y=A[j+k+mid]*W%mod;
A[j+k]=(x+y)%mod;
A[j+k+mid]=(x-y+mod)%mod;
}
}
}
}
signed main(){
inv3=powi(3,mod-2);
//cout<<inv3<<endl;
for(int i=1;i<=600;i<<=1)
powi3[i]=powi(3,(mod-1)/i),powinv3[i]=powi(inv3,(mod-1)/i);
while(getline(cin,t)){
limit=1;len=0;
memset(rev,0,sizeof(rev));
memset(a,0,sizeof(a));
memset(b,0,sizeof(b));
int tot=0,now=0,cheng=0;
int lena=0,lenb=0;
//printf("%s\n",t);
for(int i=0;i<(int)t.length();){
//cout<<t[i]<<endl;
if(t[i]=='('){
//cout<<"Case 1 start "<<endl;
memset(s,0,sizeof(s));
tot=0;++now;
i++;
}else if(t[i]==')'){
//cout<<"Case 2 endi "<<endl;
sort(s+1,s+1+tot);
//cout<<"tot "<<tot<<endl;
//for(int j=1;j<=tot;j++)cout<<s[j].x<<" "<<s[j].y<<endl;
for(int j=1;j<=tot;j++)
if(now==1)a[s[j].y]=s[j].x,lena=max(lena,s[j].y);
else if(now==2)b[s[j].y]=s[j].x,lenb=max(lenb,s[j].y);
i++;
}else if(t[i]=='*')i++,++cheng;
else if(t[i]==' ')i++;
else{
//cout<<"Case 3 mid "<<endl;
int x=0;
while(t[i]>='0'&&t[i]<='9'){
x=x*10+t[i]-'0';i++;
}
++tot;
s[tot].x=x;
if(t[i]!=')'){
i+=2;
x=0;
while(t[i]>='0'&&t[i]<='9'){
x=x*10+t[i]-'0';i++;
}
s[tot].y=x;
}
if(t[i]=='+')i++;
}
}
if(!cheng){
cout<<t<<endl;
continue;
}
//cout<<lena<<" "<<lenb<<endl;
//cout<<"lena "<<lena<<" lenb "<<lenb<<endl;
//cout<<"a "<<endl;
//for(int i=0;i<=lena;i++)cout<<a[i]<<" ";cout<<endl;
//cout<<"b "<<endl;
//for(int i=0;i<=lenb;i++)cout<<b[i]<<" ";cout<<endl;
while(limit<=lena+lenb)limit<<=1,++len;
for(int i=0;i<limit;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
//cout<<"Case 1 "<<endl;
NTT(a,1);NTT(b,1);
//cout<<"Case 2 "<<endl;
for(int i=0;i<limit;i++)a[i]=a[i]*b[i]%mod;
NTT(a,-1);
inv=powi(limit,mod-2);
for(int i=0;i<=lena+lenb;i++)a[i]=a[i]*inv%mod;
//for(int i=0;i<=lena+lenb;i++)cout<<a[i]<<" ";cout<<endl;
int n=lena+lenb+2;
while(a[n]==0)--n;
cheng=0;
while(a[cheng]==0)++cheng;
for(int i=n;i>=0;i--)
if(i>0&&a[i]>0){
printf("%da^%d",a[i],i);
if(cheng!=i)printf("+");
}else if(i==0&&a[i])printf("%d",a[i]);
cout<<endl;
}
}
P3338 [ZJOI2014]力（FTT模板1.0）
https://www.luogu.com.cn/problem/P3338
不知道为什么，重构一遍就过了……
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<iostream>
#include<cmath>
#include<bits/stdc++.h>
using namespace std;
const int N=4e5+10;
const double pi=acos(-1.0);
int n,limit=1,len,rev[N];
struct node{
double x,y;
node(double xi=0,double yi=0){
x=xi;y=yi;
}
node operator +(const node &b)const{
return {x+b.x,y+b.y};
}
node operator -(const node &b)const{
return {x-b.x,y-b.y};
}
node operator *(const node &b)const{
return {x*b.x-y*b.y,x*b.y+y*b.x};
}
}H[N],G[N],I[N];
inline void FFT(node *A,double flag){
for(int i=0;i<limit;i++)if(i<rev[i])swap(A[i],A[rev[i]]);
for(int mid=1;mid<limit;mid<<=1){
node Wn={cos(pi/mid),flag*sin(pi/mid)};
for(int j=0;j<limit;j+=(mid<<1)){
node W={1.0,0};
for(int k=0;k<mid;k++,W=W*Wn){
node x=A[j+k],y=A[j+k+mid]*W;
A[j+k]=x+y;
A[j+k+mid]=x-y;
}
}
}
}
signed main(){
cin>>n;
for(int i=1;i<=n;i++){
cin>>H[i].x;
G[n-i+1].x=H[i].x;
I[i].x=(double)1/(double)i/(double)i;
}
while(limit<=(n<<1))limit<<=1,++len;
for(int i=0;i<limit;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
FFT(H,1);FFT(G,1);FFT(I,1);
for(int i=0;i<=limit;i++)H[i]=H[i]*I[i],G[i]=G[i]*I[i];
FFT(H,-1);FFT(G,-1);
for(int i=0;i<=limit;i++)H[i].x/=limit,G[i].x/=limit;
for(int i=1;i<=n;i++)printf("%.4lf\n",H[i].x-G[n-i+1].x);
return 0;
}
\[\lfloor \quad \rfloor\\
\]
posted on
浙公网安备 33010602011771号