题目描述 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了
1 for i:=0 to n>>(t+1)-1 do 2 begin 3 p:=i<<(t+1)+s; 4 wt:=w[i<<t]*a[p+1<<t]; 5 tt[i]:=a[p]+wt; 6 tt[i+n>>(t+1)]:=a[p]-wt; 7 end;
先把单位根的几个性质搞懂,然后就基本上可以理解这个分治了(最好是自己手算一下......我就懒得说什么了,自己都不是很理解,代码还是抄别人的)
1 const 2 s=4; 3 type 4 cp=record 5 x,y:double; 6 end; 7 arr=array[0..1 shl 17]of cp; 8 var 9 c:array[0..100000]of int64; 10 a,b,w,tt:arr; 11 n,i,p:longint; 12 st:ansistring; 13 wt:cp; 14 15 operator *(var a,b:cp)c:cp; 16 begin 17 c.x:=a.x*b.x-a.y*b.y; 18 c.y:=a.x*b.y+a.y*b.x; 19 end; 20 21 operator +(var a,b:cp)c:cp; 22 begin 23 c.x:=a.x+b.x; 24 c.y:=a.y+b.y; 25 end; 26 27 operator -(var a,b:cp)c:cp; 28 begin 29 c.x:=a.x-b.x; 30 c.y:=a.y-b.y; 31 end; 32 33 procedure dft(var a:arr;s,t:longint); 34 begin 35 if n>>t=1 then exit; 36 dft(a,s,t+1); 37 dft(a,s+1<<t,t+1); 38 for i:=0 to n>>(t+1)-1 do 39 begin 40 p:=i<<(t+1)+s; 41 wt:=w[i<<t]*a[p+1<<t]; 42 tt[i]:=a[p]+wt; 43 tt[i+n>>(t+1)]:=a[p]-wt; 44 end; 45 for i:=0 to n>>t-1 do 46 a[i<<t+s]:=tt[i]; 47 end; 48 49 procedure get(var a:arr); 50 var 51 i,l,ll:longint; 52 c:char; 53 begin 54 read(c); 55 st:=''; 56 while (ord(c)>=ord('0')) and (ord(c)<=ord('9')) do 57 begin 58 st:=st+c; 59 read(c); 60 end; 61 while length(st)mod s<>0 do 62 insert('0',st,0); 63 ll:=length(st)div s; 64 for l:=1 to ll do 65 val(copy(st,l*s-s+1,s),a[ll-l].x); 66 while n>>1<l do 67 n:=n<<1; 68 end; 69 70 begin 71 n:=1; 72 get(a); 73 get(b); 74 for i:=0 to n-1 do 75 w[i].x:=cos(pi*2*i/n); 76 for i:=0 to n-1 do 77 w[i].y:=sin(pi*2*i/n); 78 dft(a,0,0); 79 dft(b,0,0); 80 for i:=0 to n-1 do 81 a[i]:=a[i]*b[i]; 82 for i:=0 to n-1 do 83 w[i].y:=-w[i].y; 84 dft(a,0,0); 85 for i:=0 to n-1 do 86 begin 87 c[i]:=c[i]+round(a[i].x/n); 88 c[i+1]:=c[i] div 10000; 89 c[i]:=c[i] mod 10000; 90 end; 91 while (c[i]=0) and (i>0) do 92 dec(i); 93 for p:=i downto 0 do 94 begin 95 str(c[p],st); 96 while (i<>p) and (length(st)<s) do 97 st:='0'+st; 98 write(st); 99 end; 100 end.
又写了一遍233
1 const 2 maxn=400400; 3 type 4 cp=record 5 x,y:double; 6 end; 7 arr=array[0..maxn]of cp; 8 var 9 a,b,w,tt:arr; 10 n:longint; 11 s:ansistring; 12 c:array[0..maxn]of longint; 13 14 operator -(a,b:cp)c:cp; 15 begin 16 c.x:=a.x-b.x; 17 c.y:=a.y-b.y; 18 end; 19 20 operator +(a,b:cp)c:cp; 21 begin 22 c.x:=a.x+b.x; 23 c.y:=a.y+b.y; 24 end; 25 26 operator *(a,b:cp)c:cp; 27 begin 28 c.x:=a.x*b.x-a.y*b.y; 29 c.y:=a.x*b.y+a.y*b.x; 30 end; 31 32 procedure get(var a:arr); 33 var 34 c:char; 35 i,len:longint; 36 begin 37 s:=''; 38 read(c); 39 while c in['0'..'9'] do 40 begin 41 s:=s+c; 42 read(c); 43 end; 44 len:=length(s); 45 for i:=1 to len do 46 a[len-i].x:=ord(s[i])-ord('0'); 47 while n<len<<1 do 48 n:=n<<1; 49 end; 50 51 procedure dft(var a:arr;s,t:longint); 52 var 53 i,p:longint; 54 wt:cp; 55 begin 56 if n>>t=1 then exit; 57 dft(a,s+1<<t,t+1);dft(a,s,t+1); 58 for i:=0 to n>>t>>1-1 do 59 begin 60 p:=i<<t<<1+s; 61 wt:=w[i<<t]*a[p+1<<t]; 62 tt[i]:=a[p]+wt; 63 tt[i+n>>t>>1]:=a[p]-wt; 64 end; 65 for i:=0 to n>>t-1 do 66 a[s+i<<t]:=tt[i]; 67 end; 68 69 procedure main; 70 var 71 i:longint; 72 begin 73 n:=1;get(a);get(b); 74 for i:=0 to n-1 do w[i].x:=cos(pi*2*i/n); 75 for i:=0 to n-1 do w[i].y:=sin(pi*2*i/n); 76 dft(a,0,0);dft(b,0,0); 77 for i:=0 to n-1 do a[i]:=a[i]*b[i]; 78 for i:=0 to n-1 do w[i].y:=-w[i].y; 79 dft(a,0,0); 80 for i:=0 to n-1 do c[i]:=round(a[i].x/n); 81 for i:=0 to n-2 do 82 begin 83 inc(c[i+1],c[i]div 10); 84 c[i]:=c[i]mod 10; 85 end; 86 i:=n-1; 87 while (c[i]=0) and (i>0) do dec(i); 88 for i:=i downto 0 do write(c[i]); 89 end; 90 91 begin 92 main; 93 end.