大数乘法
1 /* 2 题意:大数乘法NTT 3 题解:FNT 4 时间:2018.07.30 5 */ 6 7 #include <bits/stdc++.h> 8 using namespace std; 9 10 typedef long long LL; 11 const int MAXN = 100005; 12 const LL MOD7 = 1e9+7; 13 const LL MOD717 = 998244353LL; // 7*17*2^23 +1 14 const LL MOD479 = (479LL<<21) + 1; // 1004535809 15 16 inline LL Add(LL a,LL b, LL P) { a+=b;if(a>=P) a-=P;return a; } 17 inline LL Del(LL a,LL b, LL P) { a-=b;if(a<0LL) a+=P;return a; } 18 inline LL Mul(LL a,LL b, LL P) { a=a*b%P;return a; } 19 20 LL Pow(LL a,LL b, LL P) 21 { 22 LL res=1LL; 23 LL ans=a%P; 24 while (b) 25 { 26 if (b&1) res=res*ans%P; 27 ans=ans*ans%P; 28 b>>=1; 29 } 30 return res; 31 } 32 33 34 struct NTT 35 { 36 int L,l,rev[MAXN*20]; 37 38 int pre(int n) 39 { 40 for (l=0;(1<<l)<n;++l); 41 L = (1<<l); 42 for (int i=0;i<L;++i) 43 rev[i] = rev[i >> 1] >> 1 | (i&1) << (l-1); 44 return L; 45 } 46 47 void dft(LL *a, LL P) 48 { 49 //蝴蝶操作交换每个系数到需要的位置上 50 for (int i=0;i<L;++i) 51 { 52 if (i>rev[i]) swap(a[i],a[rev[i]]); 53 } 54 // 从下向上计算,i枚举宽度,第一层不需要计算,从步长为2的开始 55 // m 是步长 为宽度的一半 56 // wn 是第n层的单位根 (wn)^n %P = 1, 3是模P的原根, L | (P-1) 保证原根满足FFT的单位根的性质 57 // j 枚举每个宽度的开始位置 58 // k 枚举每个宽度内的位置 59 // w 为每次蝴蝶操作的系数 60 // 每次操作的两个位置为a[j+k] 和a[j+k+m] 并且a[j+k] = a[j+k] + w*a[j+k+m] a[j+k+m] = a[j+k] - w*a[j+k+m] 61 for (int i=2;i<=L;i<<=1) 62 { 63 int m=(i>>1); 64 LL wn = Pow(3, (P-1)/i, P); 65 // printf("wn=%lld\n",wn); 66 for (int j=0;j<L;j+=i) 67 { 68 LL w=1LL; 69 for (int k=0;k!=m;++k) 70 { 71 LL u = a[j+k], v=Mul(a[j+k+m],w,P); 72 // printf("m+j+k=%d u=%lld v=%lld\n",m+j+k,u,v); 73 a[j+k+m] = Del(u,v,P); 74 a[j+k] = Add(u,v,P); 75 w=Mul(w, wn, P); 76 } 77 } 78 } 79 } 80 81 void idft(LL *a,LL P) 82 { 83 // FFT 的逆变换需要-wn来控制, 这里通过reverse操作来替换 84 // xn = 1/N sum(Xk * g^(-k)) 85 // = 1/N sum(Xk * g^(-k) * g^n) g^n % P = 1 86 // = 1/N sum(XK * g^(n-k)) 87 // 而 X0 的系数变成 X0 * g^n = X0 保持不动, 而后面的位置需要逆序 88 reverse(a+1,a+L); 89 dft(a,P); 90 LL inv=Pow(L,P-2,P); 91 for (int i=0;i<L;++i) a[i]=a[i]*inv%P; 92 } 93 }ntt; 94 95 void Print(vector<LL> a) 96 { 97 for (auto it=a.begin();it!=a.end();++it) 98 { 99 printf("%lld ",*it); 100 } 101 puts(""); 102 } 103 104 105 vector<LL> Mul(vector<LL> a, vector<LL> b, LL P) 106 { 107 int n = a.size() + b.size() - 1; 108 // printf("n=%d\n",n); 109 int L = ntt.pre(n); 110 // printf("L=%d\n",L); 111 vector<LL> c; 112 // printf("L=%d\n",L); 113 114 if (1LL * a.size()*b.size() <= (a.size() + b.size())*20LL) 115 { 116 c.resize(n); 117 for (int i=0;i<a.size();++i) 118 { 119 for (int j=0;j<b.size();++j) 120 { 121 c[i+j] += a[i] * b[j] % P; 122 } 123 } 124 return c; 125 } 126 c.resize(L);a.resize(L),b.resize(L); 127 ntt.dft(&a[0],P); 128 // printf("ntt a\n"); 129 // Print(a); 130 ntt.dft(&b[0],P); 131 //////printf("ntt b\n"); 132 // Print(b); 133 for (int i=0;i<L;++i) 134 c[i] = a[i]*b[i]%P; 135 ntt.idft(&c[0],P); 136 c.resize(n); 137 // printf("idft c\n"); 138 // Print(c); 139 return c; 140 } 141 142 char sa[MAXN]; 143 char sb[MAXN]; 144 145 vector<LL> a,b,c; 146 vector<LL> ans; 147 148 149 150 151 int main() 152 { 153 #ifndef ONLINE_JUDGE 154 freopen("test.txt","r",stdin); 155 #endif // ONLINE_JUDGE 156 while (scanf("%s%s",sa,sb)!=-1) 157 { 158 a.clear();b.clear();c.clear();ans.clear(); 159 int i,j,n; 160 char *s; 161 s=sa; 162 n = strlen(s); 163 for (j=0;j<n && s[j]=='0';++j); 164 if (j==n) --j; 165 for (i=n-1;i>=j;--i) a.push_back(s[i]-'0'); 166 s = sb; 167 n=strlen(s); 168 for (j=0;j<n&&s[j]=='0';++j); 169 if (j==n) --j; 170 for (i=n-1;i>=j;--i) b.push_back(s[i]-'0'); 171 c = Mul(a,b,MOD717); 172 LL tmp,ci=0; 173 for (i=0;i<c.size();++i) 174 { 175 tmp=c[i]+ci; 176 ans.push_back(tmp%10); 177 ci=tmp/10; 178 } 179 while (ci) 180 { 181 ans.push_back(ci%10); 182 ci/=10; 183 } 184 i=ans.size()-1; 185 for (;i>0 && ans[i]==0;--i,ans.pop_back()); 186 reverse(ans.begin(),ans.end()); 187 for (auto it=ans.begin();it!=ans.end();++it) 188 { 189 printf("%lld",*it); 190 } 191 printf("\n"); 192 } 193 return 0; 194 }