3123 高精度练习之超大整数乘法 - Wikioi

题目描述 Description

给出两个正整数A和B，计算A*B的值。保证A和B的位数不超过100000位。

输入描述 Input Description

读入两个用空格隔开的正整数

输出描述 Output Description

输出A*B的值

样例输入 Sample Input

4 9

样例输出 Sample Output

36

数据范围及提示 Data Size & Hint

两个正整数的位数不超过100000位

什么都不说了，裸FFT

看了一天的算导，勉强看懂了怎么O(nlogn)求dft，求dft^(-1)就实在看不懂了，唉，只好记下结论

首先是理解系数表达和点值表达，因为点值表达乘法是O（n）的，所以我们要通过系数表达算出点值表达，相乘后再算出点值表达

这个是分治求的dft奇偶分成两个序列，分别求出这两个序列的dft然后通过下面的代码就可以O（n）合并这个两个dft了

for i:=0 to n>>(t+1)-1 do
begin
p:=i<<(t+1)+s;
wt:=w[i<<t]*a[p+1<<t];
tt[i]:=a[p]+wt;
tt[i+n>>(t+1)]:=a[p]-wt;
end;

先把单位根的几个性质搞懂，然后就基本上可以理解这个分治了（最好是自己手算一下......我就懒得说什么了，自己都不是很理解，代码还是抄别人的）

zky神犇手动搬运算导FFT

pascal代码（我从这里抄的）

const
maxn=400400;
type
cp=record
x,y:double;
end;
arr=array[0..maxn]of cp;
var
a,b,w,tt:arr;
n:longint;
s:ansistring;
c:array[0..maxn]of longint;

operator -(a,b:cp)c:cp;
begin
c.x:=a.x-b.x;
c.y:=a.y-b.y;
end;

operator +(a,b:cp)c:cp;
begin
c.x:=a.x+b.x;
c.y:=a.y+b.y;
end;

operator *(a,b:cp)c:cp;
begin
c.x:=a.x*b.x-a.y*b.y;
c.y:=a.x*b.y+a.y*b.x;
end;

procedure get(var a:arr);
var
c:char;
i,len:longint;
begin
s:='';
read(c);
while c in['0'..'9'] do
begin
s:=s+c;
read(c);
end;
len:=length(s);
for i:=1 to len do
a[len-i].x:=ord(s[i])-ord('0');
while n<len<<1 do
n:=n<<1;
end;

procedure dft(var a:arr;s,t:longint);
var
i,p:longint;
wt:cp;
begin
if n>>t=1 then exit;
dft(a,s+1<<t,t+1);dft(a,s,t+1);
for i:=0 to n>>t>>1-1 do
begin
p:=i<<t<<1+s;
wt:=w[i<<t]*a[p+1<<t];
tt[i]:=a[p]+wt;
tt[i+n>>t>>1]:=a[p]-wt;
end;
for i:=0 to n>>t-1 do
a[s+i<<t]:=tt[i];
end;

procedure main;
var
i:longint;
begin
n:=1;get(a);get(b);
for i:=0 to n-1 do w[i].x:=cos(pi*2*i/n);
for i:=0 to n-1 do w[i].y:=sin(pi*2*i/n);
dft(a,0,0);dft(b,0,0);
for i:=0 to n-1 do a[i]:=a[i]*b[i];
for i:=0 to n-1 do w[i].y:=-w[i].y;
dft(a,0,0);
for i:=0 to n-1 do c[i]:=round(a[i].x/n);
for i:=0 to n-2 do
begin
inc(c[i+1],c[i]div 10);
c[i]:=c[i]mod 10;
end;
i:=n-1;
while (c[i]=0) and (i>0) do dec(i);
for i:=i downto 0 do write(c[i]);
end;

begin
main;
end.

View Code