CODEVS3123 a*b problem plus (FFT)
1 type xh=record 2 x,y:double; 3 end; 4 arr=array[0..1000008] of xh; 5 var n,m:longint; 6 s1,s2:ansistring; 7 a,b,g,w:arr; 8 ch:char; 9 operator -(a,b:xh) c:xh; 10 begin 11 c.x:=a.x-b.x; 12 c.y:=a.y-b.y; 13 end; 14 operator +(a,b:xh) c:xh; 15 begin 16 c.x:=a.x+b.x; 17 c.y:=a.y+b.y; 18 end; 19 operator *(a,b:xh) c:xh; 20 begin 21 c.x:=a.x*b.x-a.y*b.y; 22 c.y:=a.x*b.y+a.y*b.x; 23 end; 24 procedure dft(var a:arr;s,t:longint);//a待处理数组 s初始位置 1<<t长度 25 var i,p:longint; 26 cnt:xh; 27 begin 28 if n>>t=1 then exit; 29 dft(a,s,t+1); dft(a,s+1<<t,t+1); 30 for i:=0 to n>>t>>1-1 do 31 begin 32 p:=i<<t<<1+s; 33 cnt:=w[i<<t]*a[p+1<<t]; 34 g[i]:=a[p]+cnt; 35 g[i+n>>t>>1]:=a[p]-cnt; 36 end; 37 for i:=0 to n>>t-1 do a[s+i<<t]:=g[i]; 38 end; 39 procedure clr(var a:arr); 40 begin 41 fillchar(a,sizeof(a),0); 42 end; 43 procedure FFT; 44 var i,len:longint; 45 lx:longint; 46 begin 47 n:=1; 48 while n<m<<1 do n:=n<<1; 49 for i:=0 to n-1 do w[i].x:=cos(pi*2*i/n); 50 for i:=0 to n-1 do w[i].y:=sin(pi*2*i/n); 51 dft(a,0,0); dft(b,0,0); 52 for i:=0 to n-1 do a[i]:=a[i]*b[i]; 53 for i:=0 to n-1 do w[i].y:=-w[i].y; //!!!! 54 dft(a,0,0); 55 for i:=m<<1-1 downto 0 do a[i].x:=a[i].x/n;//!!!! 56 for i:=0 to m<<1-1 do 57 begin 58 lx:=round(a[i].x); 59 a[i+1].x:=a[i+1].x+lx div 10; 60 a[i].x:=lx mod 10; 61 end; 62 len:=m<<1-1; 63 while a[len].x<1e-12 do dec(len); 64 for i:=len downto 0 do write(a[i].x:0:0); 65 writeln; 66 end; 67 procedure main; 68 var i:longint; 69 begin 70 read(ch); 71 while ord(ch)>=48 do begin s1:=s1+ch; read(ch); end; 72 readln(s2); 73 clr(a); clr(b); 74 for i:=length(s1) downto 1 do a[length(s1)-i].x:=ord(s1[i])-48; 75 for i:=length(s2) downto 1 do b[length(s2)-i].x:=ord(s2[i])-48; 76 if length(s1)>length(s2) then m:=length(s1) else m:=length(s2); 77 FFT; 78 end; 79 begin 80 main; 81 end.