大数乘法

  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 }

 

posted @ 2018-07-30 11:21  LeeSongt  阅读(219)  评论(0编辑  收藏  举报