FFT（快速傅立叶变换）

大佬博客

hdu 4609

题意:给定一个数组，问从其中选3个值能构成三角形的概率是多少。

思路:先求出选两个之和的情况，然后枚举选取的最长边，根据三角形的三边定理来求解。

选两个之和的情况及是，先将长度相同的统计起来，然后求这个数组的卷积，其值就是和为i的有多少个，也就是代码中的num数组。

/*
这里解释一下什么是两个数组的卷积，假如给定两个一元多项式，
他们两个多项式的系数会成为两个数组，而这两个多项式相乘，
会形成另一个一元多项式，而这个一元多项式的各项系数就是上面
两个数组的卷积。
*/
#include <cmath>//fft求a数组和b数组的卷积
#include <cstdio>//num[i]存的就是卷积的最后结果也就是次幂为i的多项式系数是多少
#include <cstring>//比如num[2]=3代表的就是x^2前面的系数是3
#include <iostream>
#include <algorithm>
using namespace std;
#define LL long long
const double pi=acos(-1.0);
const int maxn=400040;
int x[maxn];
LL num[maxn];//存储求得的卷积结果
LL sum[maxn];//前缀和数组
struct Complex//定义一个复数类及其运算
{
double x,y;//实部，虚部
Complex(double _x=0.0,double _y=0.0)
{
x=_x;
y=_y;
}
Complex operator -(const Complex &b)const
{
return Complex(x-b.x,y-b.y);
}
Complex operator +(const Complex &b)const
{
return Complex(x+b.x,y+b.y);
}
Complex operator *(const Complex &b)const
{
return Complex(x*b.x-y*b.y,x*b.y+y*b.x);
}
} x1[maxn],x2[maxn];
void Rader(Complex y[],int len)//雷德算法--倒位序
{
int i,j,k;
for(i=1,j=len/2; i<len-1; i++)
{
if(i<j) swap(y[i],y[j]);
k=len/2;
while(j>=k)
{
j-=k;
k/=2;
}
if(j<k) j+=k;
}
}
void FFT(Complex y[],int len,int on)//FFT实现
{
Rader(y,len);
for(int h=2; h<=len; h<<=1)
{
Complex wn(cos(-on*2*pi/h),sin(-on*2*pi/h));
for(int j=0; j<len; j+=h)
{
Complex w(1,0);//旋转因子
for(int k=j; k<j+h/2; k++)
{
Complex u=y[k];
Complex t=w*y[k+h/2];
y[k]=u+t;//蝴蝶合并操作
y[k+h/2]=u-t;
w=w*wn;//更新旋转因子
}
}
}
if(on==-1)
{
for(int i=0; i<len; i++)
y[i].x/=len;
}
}
void Conv(Complex a[],Complex b[],int len)//求卷积
{
FFT(a,len,1);
FFT(b,len,1);//求值 将系数表示转化为点值表示
for(int i=0; i<len; i++)
a[i]=a[i]*b[i];
FFT(a,len,-1);//插值 将点值表示转化为系数表示
for(int i=0; i<len; i++)
num[i]=(LL)(a[i].x+0.5);
}
int init(LL a[],LL b[],int n)
{
int len1=x[n-1]+1,len2=x[n-1]+1;
int len=1;//最后的总长度
while(len<2*len1||len<2*len2) len<<=1;
for(int i=0; i<len1; i++)
x1[i]=Complex(a[i],0);
for(int i=len1; i<len; i++)
x1[i]=Complex(0,0);
for(int i=0; i<len2; i++)
x2[i]=Complex(b[i],0);
for(int i=len2; i<len; i++)
x2[i]=Complex(0,0);
return len;
}
int main()
{
int ncase;
scanf("%d",&ncase);
while(ncase--)
{
memset(num,0,sizeof(num));
int n;
scanf("%d",&n);
for(int i=0; i<n; i++)
{
scanf("%d",&x[i]);
num[x[i]]++;
}
sort(x,x+n);
int len=init(num,num,n);
Conv(x1,x2,len);
for(int i=0; i<n; i++) //减去2次选的同一个数
num[x[i]+x[i]]--;
for(int i=1; i<=len; i++) num[i]/=2; //选1 2和选2 1是一样的所以除2
sum[0]=0;
for(int i=1; i<=len; i++) //num的前缀和
sum[i]=sum[i-1]+num[i];
LL ans=0;
for(int i=0; i<n; i++) //取3个a[i]为最长边的取法
{
ans+=sum[len]-sum[x[i]];//a[i]是最大的所有情况
ans-=(LL)(n-1-i)*i;//一个比a[i]大一个比a[i]小
ans-=(LL)(n-1);//一个取a[i]另一个取其他的
ans-=(LL)(n-1-i)*(n-2-i)/2;//两个都比a[i]大
}
LL all=1ll*n*(n-1)*(n-2)/6;//所有情况
printf("%.7lf\n",1.0*ans/all);
}
}

hdu 1402

题意:大数乘法卡时间要nlognFFT优化

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define LL long long
const double pi=acos(-1.0);
const int maxn=500050;
char s1[maxn],s2[maxn];
int num[maxn];//存储求得的卷积结果
struct Complex//定义一个复数类及其运算
{
double x,y;//实部，虚部
Complex(double _x=0.0,double _y=0.0)
{
x=_x;
y=_y;
}
Complex operator -(const Complex &b)const
{
return Complex(x-b.x,y-b.y);
}
Complex operator +(const Complex &b)const
{
return Complex(x+b.x,y+b.y);
}
Complex operator *(const Complex &b)const
{
return Complex(x*b.x-y*b.y,x*b.y+y*b.x);
}
} x1[maxn],x2[maxn];
void Rader(Complex y[],int len)//雷德算法--倒位序
{
int i,j,k;
for(i=1,j=len/2; i<len-1; i++)
{
if(i<j) swap(y[i],y[j]);
k=len/2;
while(j>=k)
{
j-=k;
k/=2;
}
if(j<k) j+=k;
}
}
void FFT(Complex y[],int len,int on)//FFT实现
{
Rader(y,len);
for(int h=2; h<=len; h<<=1)
{
Complex wn(cos(-on*2*pi/h),sin(-on*2*pi/h));
for(int j=0; j<len; j+=h)
{
Complex w(1,0);//旋转因子
for(int k=j; k<j+h/2; k++)
{
Complex u=y[k];
Complex t=w*y[k+h/2];
y[k]=u+t;//蝴蝶合并操作
y[k+h/2]=u-t;
w=w*wn;//更新旋转因子
}
}
}
if(on==-1)
{
for(int i=0; i<len; i++)
y[i].x/=len;
}
}
int len;
void Conv()//求卷积
{
FFT(x1,len,1);
FFT(x2,len,1);//求值 将系数表示转化为点值表示
for(int i=0; i<len; i++)
x1[i]=x1[i]*x2[i];
FFT(x1,len,-1);//插值 将点值表示转化为系数表示
for(int i=0; i<len; i++)
num[i]=x1[i].x+0.5;
}
void init()
{
int len1=strlen(s1),len2=strlen(s2);
len=1;
while(len<2*len1||len<2*len2) len<<=1;
for(int i=0; i<len1; i++)
x1[i]=Complex(s1[len1-i-1]-'0',0.0);
for(int i=len1; i<len; i++)
x1[i]=Complex(0.0,0.0);
for(int i=0; i<len2; i++)
x2[i]=Complex(s2[len2-i-1]-'0',0.0);
for(int i=len2; i<len; i++)
x2[i]=Complex(0.0,0.0);
}
void solve()
{
for(int i=0;i<=len;i++)
{
num[i+1]+=num[i]/10;
num[i]%=10;
}
int f=0;
for(int i=len;i>=0;i--)
{
if(num[i])
{
for(int j=i;j>=0;j--)
f=1,printf("%d",num[j]);
break;
}
}
if(!f) printf("0");
printf("\n");
}
int main()
{
while(~scanf("%s%s",s1,s2))
{
memset(num,0,sizeof(num));
init();
Conv();
solve();
}
}