【快速傅里叶变换】FFT高精乘法

比较优美的算法,虽然我至今还推倒不出那个定理(多项式差值的唯一性还有逆运算的公式),嘛,记住算了。

 

可以用来优化高精度乘法,做到NlogN,不过这个NlogN真的很蛋疼,我四位一压才可以做到300000一秒出,因为涉及复数和实数运算,常数太大了。

 

基本上是照着算导写的,常数小优化就是那些复数根可以预先处理出来,代码小优化就是DFT和DFT-1可以写在一起,只要把单位复根的y取反即可。

 

看了一下h8oj发现rank1的神牛比我快3倍,难道是底层优化?

rank4。。。

 

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

 

点赞