码迷,mamicode.com
首页 > 其他好文 > 详细

再写FFT模板

时间:2015-05-26 23:12:31      阅读:198      评论:0      收藏:0      [点我收藏+]

标签:

  没什么好说的,今天有考了FFT(虽然不用FFT也能过)但是确实有忘了怎么写FFT了,于是乎只有重新写一遍FFT模板练一下手了。第一部分普通FFT,第二部分数论FFT,记一下模数2^23*7*17+1

 

  

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
typedef double real;
namespace task1
{
        const int maxn = 110000;
        const real pi = 3.1415926535897932;
        struct Complex
        {
                real x,y;
                Complex(){}
                Complex(real x,real y=0):x(x),y(y){}
        };
        Complex operator + (Complex c1,Complex c2)
        {
                return Complex(c1.x+c2.x,c1.y+c2.y);
        }
        Complex operator - (Complex c1,Complex c2)
        {
                return Complex(c1.x-c2.x,c1.y-c2.y);
        }
        Complex operator * (Complex c1,Complex c2)
        {
                return Complex(c1.x*c2.x-c1.y*c2.y,c1.y*c2.x+c1.x*c2.y);
        }
        Complex operator / (Complex c1,real k)
        {
                return Complex(c1.x/k,c1.y/k);
        }
        char s1[maxn],s2[maxn];
        Complex g1[maxn],g2[maxn],rr[maxn];
        int lst[maxn];
        void FFT(Complex g[],int l,int p)
        {
                int x;
                for (int i=0;i<l;i++)
                {
                        x=0;
                        for (int j=1,k=(l>>1);j<l;j<<=1,k>>=1)
                                if (i&j)x+=k;
                        if (x<i)
                                swap(g[x],g[i]);
                }
                Complex wt,w,tmp;
                for (int i=1;i<l;i<<=1)
                {
                        wt=Complex(cos(pi/i*p),sin(pi/i*p));
                        for (int j=0;j<l;j+=(i<<1))
                        {
                                Complex w=Complex(1,0);
                                for (int k=0;k<i;k++)
                                {
                                        tmp=g[j+k];
                                        g[j+k]=g[j+k]+g[i+j+k]*w;
                                        g[i+j+k]=tmp-g[i+j+k]*w;
                                        w=w*wt;
                                }
                        }
                }
        }
        void main()
        {
                int l1,l2;
                scanf("%s",s1);
                scanf("%s",s2);
                l1=(int)strlen(s1);
                l2=(int)strlen(s2);
                for (int i=0;i<l1;i++)
                        g1[(l1-i-1)/4]=g1[(l1-i-1)/4].x*10+s1[i]-0;
                for (int i=0;i<l2;i++)
                        g2[(l2-i-1)/4]=g2[(l2-i-1)/4].x*10+s2[i]-0;
                int l=max(l1,l2)/4+1;
                while (l!=(l&(-l)))l-=l&(-l);
                l<<=2;
                FFT(g1,l,1);
                FFT(g2,l,1);
                for (int i=0;i<l;i++)rr[i]=g1[i]*g2[i];
                FFT(rr,l,-1);
                for (int i=0;i<l;i++)
                {
                        rr[i]=rr[i]/l;
                        lst[i]=(int)(rr[i].x+0.5);
                }
                for (int i=0;i<l;i++)
                {
                        lst[i+1]+=lst[i]/10000;
                        lst[i]%=10000;
                }
                while (l>1 && !lst[l-1])l--;
                printf("%d",lst[l-1]);
                for (int i=l-2;i>=0;i--)
                        printf("%04d",lst[i]);
                printf("\n");
        }
}
namespace task2
{
        typedef long long qword;
        const int maxn = 10000;
        const qword p = 998244353;
        char s1[maxn],s2[maxn];
        qword g1[maxn],g2[maxn];;
        qword rr[maxn];
        qword pow_mod(qword x,qword y)
        {
                qword ret=1;
                while (y)
                {
                        if (y&1)ret=ret*x%p;
                        x=x*x%p;
                        y>>=1;
                }
                return ret;
        }
        void FFT(qword g[],int l,int f)
        {
                int x=0;
                for (int i=1;i<l;i++)
                {
                        x=0;
                        for (int j=1,k=(l>>1);j<l;j<<=1,k>>=1)
                                if (i&j)x+=k;
                        if (i<x)
                                swap(g[i],g[x]);
                }
                qword w,wt,tmp;
                for (int i=1,t=1;i<l;i<<=1,t++)
                {
                        wt=pow_mod(3,7*17*(1<<(23-t)))%p;
                        if (f==-1)
                                //wt=-wt;
                                wt=pow_mod(wt,p-2);
                        for (int j=0;j<l;j+=(i<<1))
                        {
                                w=1;
                                for (int k=0;k<i;k++)
                                {
                                        tmp=g[j+k];
                                        g[j+k]=(g[j+k]+g[i+j+k]*w)%p;
                                        g[i+j+k]=(tmp-g[i+j+k]*w)%p;
                                        w=w*wt%p;
                                }
                        }
                }
        }
        void main()
        {
                int l1,l2;
            //    for (int i=0;i<5;i++)
            //            printf("[1/%d*pi]:%d\n",(1<<i),7*17*(1<<(23-i)));
                scanf("%s",s1);
                scanf("%s",s2);
                l1=(int)strlen(s1);
                l2=(int)strlen(s2);
                for (int i=0;i<l1;i++)
                        g1[i]=s1[l1-i-1]-0;
                for (int i=0;i<l2;i++)
                        g2[i]=s2[l2-i-1]-0;
                int l=max(l1,l2);
                while (l!=(l&(-l)))l-=l&(-l);
                l<<=2;
                FFT(g1,l,1);
                FFT(g2,l,1);
                for (int i=0;i<l;i++)rr[i]=g1[i]*g2[i]%p;
                FFT(rr,l,-1);
                qword invl=pow_mod(l,p-2);
                for (int i=0;i<l;i++)rr[i]=(rr[i]*invl%p+p)%p;
                for (int i=0;i<l;i++)
                {
                        rr[i+1]+=rr[i]/10;
                        rr[i]%=10;
                }
                while (l>1 && ! rr[l-1])l--;
                for (int i=l-1;i>=0;i--)
                        printf("%d",(int)rr[i]);
        }
}

int main()
{
        freopen("input.txt","r",stdin);
        task1::main();
        task2::main();
}

 

再写FFT模板

标签:

原文地址:http://www.cnblogs.com/mhy12345/p/4531821.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!