hdu1402A * B Problem Plus(NTT)

来源：互联网发布：查看端口是否的命令编辑：程序博客网时间：2024/06/06 07:24

链接：http://acm.hdu.edu.cn/showproblem.php?pid=1402

题意：大数乘法

分析：练习下NTT。学习网站：http://blog.csdn.net/acdreamers/article/details/39026505 http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform#i-15

代码：

#include<map>#include<set>#include<cmath>#include<queue>#include<bitset>#include<math.h>#include<cstdio>#include<vector>#include<string>#include<cstring>#include<iostream>#include<algorithm>#pragma comment(linker, "/STACK:102400000,102400000")using namespace std;const int N=50010;const int MAX=151;const int mod=100000000;const int MOD1=100000007;const int MOD2=100000009;const double EPS=0.00000001;typedef long long ll;const ll MOD=1000000009;const ll INF=10000000010;typedef double db;typedef unsigned long long ull;ll G=3,P=469762049,a[3*N],b[3*N],wn[25];char s[N],t[N];ll q_pow(ll a,ll b,ll M) {    ll ret=1;a%=M;    while (b) {        if (b&1) ret=ret*a%M;        a=a*a%M;        b>>=1;    }    return ret;}void getwn() {    for (int i=0;i<25;i++) wn[i]=q_pow(G,(P-1)/(1<<i),P);}void NTT(ll x[],int n,int rev) {    int i,j,k,t,ds,id=0;    ll w,u,v;    for (i=1,j=n>>1,k=n>>1;i<n-1;i++,k=n>>1) {//Rader        if (i<j) swap(x[i],x[j]);        while (j>=k) { j-=k;k>>=1; }        if (j<k) j+=k;    }    for (i=2,ds=1;i<=n;i<<=1,ds++)        for (j=0;j<n;j+=i) {            w=1;            for (k=j;k<j+i/2;k++) {                u=x[k]%P;v=w*x[k+i/2]%P;                x[k]=(u+v)%P;                x[k+i/2]=(u-v+P)%P;                w=w*wn[ds]%P;            }        }    if (rev==-1) {        for (i=1;i<n/2;i++) swap(x[i],x[n-i]);        w=q_pow(n,P-2,P);        for (i=0;i<n;i++) x[i]=x[i]*w%P;    }}void solve() {    int i,n=1,les=strlen(s),let=strlen(t);    while (n<les+let) n<<=1;    for (i=0;i<les;i++) a[i]=s[les-i-1]-'0';    for (i=les;i<n;i++) a[i]=0;    for (i=0;i<let;i++) b[i]=t[let-i-1]-'0';    for (i=let;i<n;i++) b[i]=0;    NTT(a,n,1);NTT(b,n,1);    for (i=0;i<n;i++) a[i]=a[i]*b[i]%P;    NTT(a,n,-1);    for (i=0;i<n-1;i++) a[i+1]+=a[i]/10,a[i]%=10;    while (a[n-1]>10) { a[n]+=a[n-1]/10;a[n-1]%=10;n++; }    while (n>1&&!a[n-1]) n--;    for (i=n-1;i>=0;i--) printf("%I64d", a[i]);    printf("\n");}int main(){    getwn();    while (scanf("%s%s", s, t)!=EOF) solve();    return 0;}/*常用费马素数P=r*2^k+1,g为P的原根P       r       k       g3       1       1       25       1       2       217      1    4       397      3    5       5193     3    6       5257     1    8       37681        15      9       1712289       3    12      1140961       5    13      365537       1    16      3786433      3    18      105767169     11    19      37340033     7    20      323068673        11    21      3104857601       25    22      3167772161       5       25      3469762049       7       26      31004535809      479     21      32013265921      15      27      312281701377      17      27      33221225473      3       30      575161927681     35      31      377309411329     9       33      7*/

0 0