(算法课)大整数相乘 |向量卷积|多项式相乘| 快速傅里叶变换FFT | 分治
D (1021) : 很大的AB
Time Limit: 1 Sec Memory Limit: 256 Mb Submitted: 6 Solved: 0
Description
如题,计算AB的值并输出。
Input
两行,分别代表A和B。
A和B的位数不超过 106 位。
Output
一行一个整数表示乘积。
Sample Input
1
2
Sample Output
2
Hint
一般的分治优化是不行的,分析复杂度之后再写代码。
https://www.luogu.com.cn/problem/P1919
http://81.71.47.149:20082/csgoj/contest/problem?cid=1105&pid=D
在二进制基础上修改为十进制 TLE
//在二进制基础上修改为十进制
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<vector>
typedef std::vector<int> BigBinary;
void check(BigBinary &x)
{
while(!x.empty() && x.back() == 0) x.pop_back();
int cur = 0;
for(int i = 0; i < x.size(); i ++) //进位
{
x[i] += cur;
cur = x[i] / 10; x[i] = x[i] % 10;
}
for(; cur; cur /= 10) x.push_back(cur % 10);
}
void Print(BigBinary &x)
{
check(x);
for(int i = x.size() - 1; i >= 0; i --) printf("%d", x[i]);
if(x.empty()) printf("0");
printf("\n");
}
BigBinary Add(const BigBinary &a, const BigBinary &b, int negFlag=1)
{
BigBinary c(a);
for(int i = 0; i < b.size(); i ++) c[i] += b[i] * negFlag;
return c;
}
BigBinary Minus(const BigBinary &a, const BigBinary &b)
{
return Add(a, b, -1); // 本算法保证了减法结果不小于0
}
BigBinary Mul(const BigBinary &a, const BigBinary &b)
{
BigBinary c;
c.resize(a.size() * b.size() + 1);
for(int i = 0; i < a.size(); i ++)
for(int j = 0; j < b.size(); j ++)
c[i + j] += a[i] * b[j];
return c;
}
void MulN2(BigBinary &a, int n_2)
{
int isize = a.size();
a.resize(isize + n_2);
for(int i = a.size() - 1, j = isize - 1; j >= 0; i --, j--)
a[i] = a[j];
for(int i = n_2 - 1; i >= 0; i --) a[i] = 0;
}
BigBinary FasterMul(const BigBinary &a, const BigBinary &b)
{
if(a.size() < 4) return Mul(a, b);
int n_2 = a.size() >> 1;
BigBinary A(a.begin() + n_2, a.end());
BigBinary B(a.begin(), a.begin() + n_2);
BigBinary C(b.begin() + n_2, b.end());
BigBinary D(b.begin(), b.begin() + n_2);
BigBinary A_C = FasterMul(A, C);
BigBinary B_D = FasterMul(B, D);
// AD+BC = (A+B)*(C+D)-AC-BD,该方式与课本不同,避免了减法出现负数
BigBinary ADpBC = Minus(Minus(FasterMul(Add(A, B), Add(C, D)), A_C), B_D);
MulN2(A_C, n_2 << 1); MulN2(ADpBC, n_2);
return Add(Add(A_C, ADpBC), B_D);
}
const int maxn = 1e6 + 10;
char buf[maxn];
BigBinary a, b;
#ifdef CPC
#include<time.h>
#endif
int main()
{
while(scanf("%s", buf) != EOF)
{
a.clear();
b.clear();
for(int i = strlen(buf) - 1; i >= 0; i --)
a.push_back(buf[i] - '0');
scanf("%s", buf);
for(int i = strlen(buf) - 1; i >= 0; i --)
b.push_back(buf[i] - '0');
if(a.size() < b.size()) a.resize(b.size(), 0);
else b.resize(a.size(), 0);
BigBinary res = FasterMul(a, b);
Print(res);
}
return 0;
}
/*
123123333333333333333333432432423
1424222222223423424234
175354987407555303403104334472459038711630618465538982
*/
参考FFT模板https://blog.51cto.com/u_15368396/5526315
//傅里叶变换模板https://blog.51cto.com/u_15368396/5526315
//修改:注意结果进位
#include <bits/stdc++.h>
using namespace std;
const int maxn = 4e6+233;
#define DB double
const DB pi = acos(-1);
struct CP {
DB x, y; CP(){} inline CP(DB a, DB b):x(a),y(b){}
inline CP operator + (const CP&r) const { return CP(x + r.x, y + r.y); }
inline CP operator - (const CP&r) const { return CP(x - r.x, y - r.y); }
inline CP operator * (const CP&r) const { return CP(x*r.x-y*r.y, x*r.y+y*r.x); }
inline CP conj() { return CP(x, -y); }
} a[maxn], b[maxn], t;
int n, m;
inline void Swap(CP&a, CP&b) { t = a; a = b; b = t; }
void FFT(CP*a, int n, int f) {
int i, j, k;
for(i = j = 0; i < n; ++ i) {
if(i > j) Swap(a[i], a[j]);
for(k = n>>1; (j ^= k) < k; k >>= 1);
}
for(i = 1; i < n; i <<= 1) {
CP wn(cos(pi/i), f*sin(pi/i));
for(j = 0; j < n; j += i<<1) {
CP w(1, 0);
for(k = 0; k < i; ++ k, w=w*wn) {
CP x = a[j+k], y = w*a[i+j+k];
a[j+k] = x+y; a[i+j+k] = x-y;
}
}
}
if(-1 == f) for(i = 0; i < n; ++ i) a[i].x /= n;
}
char buf1[maxn];
char buf2[maxn];
int ans[maxn];
int main() {
scanf("%s", buf1);//strlen(buf)
scanf("%s", buf2);
n = strlen(buf1)-1; m = strlen(buf2)-1;
//printf("%d %d\n",n,m);
for(int i = 0; i <= n; ++ i) a[i].x = buf1[i]-'0';
for(int i = 0; i <= m; ++ i) a[i].y = buf2[i]-'0';
for(m += n, n = 1; n <= m; n <<= 1);
FFT(a, n, 1); CP Q(0, -0.25);
for(int i = 0, j; i < n; ++ i) j = (n-i)&(n-1), b[i] = (a[i]*a[i]-(a[j]*a[j]).conj())*Q;
FFT(b, n,-1);
//对结果b[i].x进行进位
//先取整
for(int i=m;i>=0;i--) ans[i]=int(b[i].x+0.2);
for(int i=m;i>=1;i--)
{
int temp=ans[i];
if(temp>9) {
ans[i]=temp%10;
ans[i-1]+=(temp/10);
}
}
for(int i = 0; i <= m; ++ i) printf("%d",ans[i]);
//fwrite(Out, 1, ou-Out, stdout);
return 0;
}
/*
123123333333333333333333432432423
1424222222223423424234
175354987407555303403104334472459038711630618465538982
*/
/**********************************************************************
Problem: 1021
User: 2210413010
Language: C++
Result: AC
Time:424 ms
Memory:150636 kb
**********************************************************************/
自己实现FFT算法 TLE
//将大整数看成向量,那么大整数相乘就是向量卷积
//卷积计算
#include <bits/stdc++.h>
using namespace std;
const int maxn = 1e6 + 10;
int ans[2*maxn+9];
char buf1[maxn];
char buf2[maxn];
int n;//向量长度
//int a[maxn]; //长度为n的向量ab,进行多项式卷积
//int b[maxn];
//选择2n个x值
const float Pi = acos(-1.0) ;
struct node{
float x, y ;
node (float xx = 0, float yy = 0){
x = xx, y = yy ;
}
}a[maxn*2+4],b[maxn*2+4],c[maxn*2+4],omiga[maxn*2+4],A[maxn*2+4],B[maxn*2+4],C[maxn*2+4],D[maxn*2+4];
node operator * (node J, node Q){
return node(J.x * Q.x - J.y * Q.y , J.x * Q.y + J.y * Q.x);
}
node operator + (node J, node Q){
return node(J.x + Q.x , J.y + Q.y);
}
node operator - (node J, node Q){
return node(J.x - Q.x , J.y - Q.y );
}
node cum1(node *F,node x,int len) //多项式求值,求F(x)
{
/**
//先用O(n)算法测试 TLE
node ans=F[len-1];
for(int i=len-2;i>=0;--i)
{
ans=x*ans+F[i];
}
return ans;
**/
//分治法多项式求值 O(logn) TLE
//递归出口
// if(len==0) return node(0,0);
if(len==1) return F[0];
//蝴蝶变换
node even[len/2+4],odd[len/2+4];
int oddi=0,eveni=0;
for(int i=1;i<len;i+=2) odd[oddi++]=F[i];
for(int i=0;i<len;i+=2) even[eveni++]=F[i];
//奇数偶数别搞反了
return cum1(even,x*x,eveni)+x*cum1(odd,x*x,oddi);
}
//https://www.cnblogs.com/lfri/p/11572792.html
node cum(node *a,int st, int en, int siz, node x)
{
// for(int i = st;i < en;i += (1 << siz)) printf("%llf ", a[i].x);
// printf("\n");
int len = (1<<siz) - 1; //间隔长度
if(st+len >= en && st <= en) return a[st];
return cum(a, st, en-len, siz+1, x*x) + x*cum(a, st+len+1, en, siz+1, x*x);
}
int main()
{
scanf("%s", buf1);
scanf("%s", buf2);//strlen(buf)
//构造向量ab
n=max(strlen(buf1),strlen(buf2));
for(int i=0;i<n;++i) {
if(((int)strlen(buf1)-i-1) >= 0)
a[n-i-1].x = (buf1[strlen(buf1)-i-1] - '0');
else
a[n-i-1].x=0;
if(((int)strlen(buf2)-i-1) >= 0)
b[n-i-1].x=(buf2[strlen(buf2)-i-1]-'0');
else
b[n-i-1].x=0;
}
//test Yes 如计算12*123 那么a 0 1 2 b 1 2 3
//确定omiga
for(int i=0;i<2*n;++i)
{
omiga[i]=node(cos(Pi*i*1.0/n),sin(Pi*i*1.0/n));
}
//FFT算法
//计算A、B、C ABC都是2n个数
for(int i=0;i<2*n;++i)
{
//A[i]=cum(a,0,n-1,0,omiga[i]);//多项式求值
//B[i]=cum(b,0,n-1,0,omiga[i]);
A[i]=cum1(a,omiga[i],n);//多项式求值
B[i]=cum1(b,omiga[i],n);
C[i]=A[i]*B[i];
}
//计算D,ci系数
for(int i=0;i<2*n;++i) {
//D[i]=cum(C,0,2*n-1,0,omiga[i]);
D[i]=cum1(C,omiga[i],2*n);
}
c[0].x=D[0].x/float(2*n);
for(int i=1;i<2*n;++i) c[2*n-i].x=D[i].x/float(2*n);
//注意进位
for(int i=0;i<2*n;++i)
ans[i]=int(c[i].x+0.2);
int ii=2*n-1;
while(1) {
if(ans[ii]==0) ii--;
else break;
}
for(int i=ii;i>=1;--i)
{
int temp=ans[i];
if(temp>9) {
ans[i]=temp%10;
ans[i-1]+=(temp/10);
}
}
//去掉前先导0
int flag1=0;
for(int i = 0; i <= ii; ++ i) {
if(ans[i]!=0) flag1=1;
if(flag1) printf("%d",ans[i]);
}
if(flag1==0) printf("0");
return 0;
}
/*
123123333333333333333333432432423
1424222222223423424234
175354987407555303403104334472459038711630618465538982
*/