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();
    }
}