您的位置:首页 > 其它

【UOJ34】高精度乘法(FFT)

2017-04-19 17:05 211 查看

题意:

思路:FFT模板,自带10倍常数

1 type cp=record
2          x,y:double;
3         end;
4      arr=array[0..510000]of cp;
5 var a,b,cur:arr;
6     n,m,n1,n2,i,j:longint;
7     x:double;
8
9 function jia(a,b:cp;f:longint):cp;
10 begin
11  if f=-1 then
12  begin
13   b.x:=-b.x; b.y:=-b.y;
14  end;
15  jia.x:=a.x+b.x;
16  jia.y:=a.y+b.y;
17 end;
18
19 function mult(a,b:cp):cp;
20 begin
21  mult.x:=a.x*b.x-a.y*b.y;
22  mult.y:=a.x*b.y+a.y*b.x;
23 end;
24
25 procedure swap(var x,y:cp);
26 var t:cp;
27 begin
28  t:=x; x:=y; y:=t;
29 end;
30
31 function max(x,y:longint):longint;
32 begin
33  if x>y then exit(x);
34  exit(y);
35 end;
36
37 procedure fft(var a:arr;n,f:longint);
38 var i,j,k,m:longint;
39     w,u,v:cp;
40
41 begin
42  i:=n>>1; j:=1;
43  while j<n do
44  begin
45   if i<j then swap(a[i],a[j]);
46   k:=n>>1;
47   while k and i>0 do
48   begin
49    i:=i xor k;
50    k:=k>>1;
51   end;
52   i:=i xor k;
53   inc(j);
54  end;
55  m:=2;
56  while m<=n do
57  begin
58   w.x:=cos(2*pi*f/m); w.y:=sin(2*pi*f/m);
59   cur[0].x:=1; cur[0].y:=0;
60   for i:=1 to m-1 do cur[i]:=mult(cur[i-1],w);
61   i:=0;
62   while i<n do
63   begin
64    j:=i;
65    while j<i+(m>>1) do
66    begin
67     u:=a[j]; v:=mult(a[j+(m>>1)],cur[j-i]);
68     a[j]:=jia(u,v,1);
69     a[j+(m>>1)]:=jia(u,v,-1);
70     inc(j);
71    end;
72    i:=i+m;
73   end;
74   m:=m<<1;
75  end;
76 end;
77
78 begin
79  assign(input,'uoj34.in'); reset(input);
80  assign(output,'uoj34.out'); rewrite(output);
81  readln(n1,n2);
82  inc(n1); inc(n2);
83  for i:=0 to n1-1 do
84  begin
85   read(x); a[i].x:=x;
86  end;
87  for i:=0 to n2-1 do
88  begin
89   read(x); b[i].x:=x;
90  end;
91  n:=max(n1,n2);
92  m:=1;
93  while m<=(n<<1) do m:=m<<1;
94  fft(a,m,1); fft(b,m,1);
95  for i:=0 to m do a[i]:=mult(a[i],b[i]);
96  fft(a,m,-1);
97  for i:=0 to n1+n2-2 do write(round(a[i].x/m),' ');
98  close(input);
99  close(output);
100 end.

 

内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: