【快速傅裡葉變換】FFT高精乘法

比較優美的演算法，雖然我至今還推倒不出那個定理（多項式差值的唯一性還有逆運算的公式），嘛，記住算了。

 

可以用來最佳化高精度乘法，做到NlogN,不過這個NlogN真的很蛋疼，我四位一壓才可以做到300000一秒出，因為涉及複數和實數運算，常數太大了。

 

基本上是照著算導寫的，常數小最佳化就是那些複數根可以預先處理出來，代碼小最佳化就是DFT和DFT-1可以寫在一起，只要把單位複根的y取反即可。

 

看了一下h8oj發現rank1的神牛比我快3倍，難道是底層最佳化？

rank4。。。

 

program fft;{$inline+}<br />const s=4;<br />type<br /> cp=record x,y:real end;<br /> arr=array[0..1 shl 17] of cp;<br />var<br /> c:array[0..100000] of int64;<br /> a,b,w,tt:arr; n,i,p:longint;<br /> st:ansistring; wt:cp;<br />operator *(var a,b:cp)c:cp;inline;<br />begin c.x:=a.x*b.x-a.y*b.y;c.y:=a.x*b.y+a.y*b.x end;<br />operator +(var a,b:cp)c:cp;inline;begin c.x:=a.x+b.x;c.y:=a.y+b.y end;<br />operator -(var a,b:cp)c:cp;inline;begin c.x:=a.x-b.x;c.y:=a.y-b.y end;<br />procedure dft(var a:arr;s,t:longint);inline;<br />begin<br /> if n>>t=1 then exit;<br /> dft(a,s,t+1);<br /> dft(a,s+1<<t,t+1);<br /> for i:=0 to n>>(t+1)-1 do begin<br /> p:=i<<(t+1)+s;<br /> wt:=w[i<<t]*a[p+1<<t];<br /> tt[i]:=a[p]+wt;<br /> tt[i+n>>(t+1)]:=a[p]-wt;<br /> end;<br /> for i:=0 to n>>t-1 do a[i<<t+s]:=tt[i];<br />end;<br />procedure get(var a:arr);var i,l,ll:longint;<br />begin<br /> readln(st);<br /> while length(st)mod s<>0 do insert('0',st,0);<br /> ll:=length(st) div s;<br /> for l:=1 to ll do val(copy(st,l*s-s+1,s),a[ll-l].x);<br /> while n>>1<l do n:=n+n;<br />end;<br />begin<br /> readln;n:=1; get(a);get(b);<br /> for i:=0 to n-1 do w[i].x:=cos(pi*2*i/n);<br /> for i:=0 to n-1 do w[i].y:=sin(pi*2*i/n);<br /> dft(a,0,0); dft(b,0,0);<br /> for i:=0 to n-1 do a[i]:=a[i]*b[i];<br /> for i:=0 to n-1 do w[i].y:=-w[i].y;<br /> dft(a,0,0);<br /> for i:=0 to n-1 do begin<br /> c[i]:=c[i]+round(a[i].x/n);<br /> c[i+1]:=c[i] div 10000;<br /> c[i]:=c[i] mod 10000;<br /> end;<br /> while (c[i]=0)and(i>0) do dec(i);<br /> for p:=i downto 0 do begin<br /> str(c[p],st);<br /> while (i<>p)and(length(st)<s) do st:='0'+st;<br /> write(st);<br /> end;<br />end.<br />