关于fft的一点总结

好吧，其实我并没有深入运用fft，只会优化卷积

听说fft经常和生成函数结合在一起………………oi真是迅猛发展，我真是与时代脱节了……

关于fft的学习推荐直接去看算法导论，写得非常清楚

主要弄懂n次单位根的相关性质定理（消去定理，折半定理）即可，当然也可以直接背代码……

bzoj2179

模板题，fft可以优化高精度乘法

顺便说一句，pascal可以定义operator，但跑得慢

这题我跑了10s……

  1 uses math;
  2 type point=record
  3        x,y:double;
  4      end;
  5      arr=array[0..150010] of point;
  6
  7 var a,b,c:arr;
  8     d,r:array[0..150010] of longint;
  9     i,n,m,l:longint;
 10     s:ansistring;
 11
 12 operator +(a,b:point)c:point;
 13   begin
 14     c.x:=a.x+b.x;
 15     c.y:=a.y+b.y;
 16   end;
 17
 18 operator -(a,b:point)c:point;
 19   begin
 20     c.x:=a.x-b.x;
 21     c.y:=a.y-b.y;
 22   end;
 23
 24 operator *(a,b:point)c:point;
 25   begin
 26     c.x:=a.x*b.x-a.y*b.y;
 27     c.y:=a.x*b.y+a.y*b.x;
 28   end;
 29
 30 procedure fft(var a:arr; ty:longint);
 31   var i,j,s,k:longint;
 32       w,p,u,v:point;
 33   begin
 34     for i:=0 to n-1 do
 35       if i<r[i] then
 36       begin
 37         w:=a[r[i]];
 38         a[r[i]]:=a[i];
 39         a[i]:=w;
 40       end;
 41     i:=1;
 42     while i<n do
 43     begin
 44       p.x:=cos(pi/i);
 45       p.y:=ty*sin(pi/i);
 46       s:=i shl 1;
 47       j:=0;
 48       while j<n do
 49       begin
 50         w.x:=1; w.y:=0;
 51         k:=0;
 52         while k<i do
 53         begin
 54           u:=a[j+k];
 55           v:=w*a[j+k+i];
 56           a[j+k]:=u+v;
 57           a[j+k+i]:=u-v;
 58           w:=w*p;
 59           inc(k);
 60         end;
 61         inc(j,s);
 62       end;
 63       i:=i shl 1;
 64     end;
 65   end;
 66
 67 begin
 68   readln(n);
 69   readln(s);
 70   for i:=0 to n-1 do
 71     a[n-1-i].x:=ord(s[i+1])-48;
 72   readln(s);
 73   for i:=0 to n-1 do
 74     b[n-1-i].x:=ord(s[i+1])-48;
 75   dec(n);
 76   m:=2*n;
 77   n:=1;
 78   l:=0;
 79   while n<=m do
 80   begin
 81     n:=n shl 1;
 82     inc(l);
 83   end;
 84   for i:=0 to n-1 do
 85     r[i]:=(r[i shr 1] shr 1) or ((i and 1) shl (l-1));
 86   fft(a,1); fft(b,1);
 87   for i:=0 to n-1 do
 88     c[i]:=a[i]*b[i];
 89   fft(c,-1);
 90   i:=0;
 91   while i<=m do
 92   begin
 93     d[i]:=d[i]+trunc(round(c[i].x/n));
 94     if d[i]>=10 then
 95     begin
 96       d[i+1]:=d[i+1]+d[i] div 10;
 97       d[i]:=d[i] mod 10;
 98     end;
 99     if d[m+1]>0 then inc(m);
100     inc(i);
101   end;
102   for i:=m downto 0 do
103     write(d[i]);
104   writeln;
105 end.

bzoj2194

模板题，倒过来就变成卷积了

不用operator就跑的快多了

 1 type point=record
 2        x,y:double;
 3      end;
 4      arr=array[0..300010] of point;
 5
 6 var a,b,c:arr;
 7     r:array[0..300010] of longint;
 8     h,i,n,l,m:longint;
 9
10 procedure fft(var a:arr; ty:longint);
11   var i,j,k,s:longint;
12       u,v,p,w:point;
13   begin
14     for i:=0 to n-1 do
15       if i<r[i] then
16       begin
17         w:=a[r[i]];
18         a[r[i]]:=a[i];
19         a[i]:=w;
20       end;
21
22     i:=1;
23     while i<n do
24     begin
25       p.x:=cos(pi/i);
26       p.y:=ty*sin(pi/i);
27       s:=i shl 1;
28       j:=0;
29       while j<n do
30       begin
31         w.x:=1;
32         w.y:=0;
33         for k:=0 to i-1 do
34         begin
35           u:=a[j+k];
36           v.x:=w.x*a[j+k+i].x-w.y*a[j+k+i].y;
37           v.y:=w.x*a[j+k+i].y+w.y*a[j+k+i].x;
38           a[j+k].x:=u.x+v.x;
39           a[j+k].y:=u.y+v.y;
40           a[j+k+i].x:=u.x-v.x;
41           a[j+k+i].y:=u.y-v.y;
42           u:=w;
43           w.x:=u.x*p.x-u.y*p.y;
44           w.y:=u.x*p.y+u.y*p.x;
45         end;
46         inc(j,s);
47       end;
48       i:=i shl 1;
49     end;
50   end;
51
52
53 begin
54   readln(n);
55   h:=n;
56   for i:=0 to n-1 do
57     readln(a[i].x,b[n-1-i].x);
58   m:=2*n-2;
59   n:=1;
60   while n<=m do
61   begin
62     n:=n*2;
63     inc(l);
64   end;
65   for i:=0 to n-1 do
66     r[i]:=(r[i shr 1] shr 1) or ((i and 1) shl (l-1));
67   fft(a,1);
68   fft(b,1);
69   for i:=0 to n-1 do
70   begin
71     c[i].x:=a[i].x*b[i].x-a[i].y*b[i].y;
72     c[i].y:=a[i].x*b[i].y+a[i].y*b[i].x;
73   end;
74   fft(c,-1);
75   for i:=h-1 to m do
76     writeln(round(c[i].x/n));
77
78 end.

bzoj3527

给一下题面吧，， 令 Ei=Fi/qi，求Ei

带进去稍微弄弄就是卷积啊……

  1 type point=record
  2        x,y:extended;
  3      end;
  4      arr=array[0..300010] of point;
  5
  6 var a,b,c:arr;
  7     ans,d:array[0..300010] of extended;
  8     r:array[0..300010] of longint;
  9     i,n,m,l:longint;
 10
 11 procedure fft(var a:arr; ty:longint);
 12   var i,j,k,s:longint;
 13       w,p,u,v:point;
 14   begin
 15     for i:=0 to n-1 do
 16       if i<r[i] then
 17       begin
 18         w:=a[r[i]];
 19         a[r[i]]:=a[i];
 20         a[i]:=w;
 21       end;
 22
 23     i:=1;
 24     while i<n do
 25     begin
 26       p.x:=cos(pi/i);
 27       p.y:=ty*sin(pi/i);
 28       s:=i shl 1;
 29       j:=0;
 30       while j<n do
 31       begin
 32         w.x:=1;
 33         w.y:=0;
 34         for k:=0 to i-1 do
 35         begin
 36           u:=a[j+k];
 37           v.x:=w.x*a[j+k+i].x-w.y*a[j+k+i].y;
 38           v.y:=w.x*a[j+k+i].y+w.y*a[j+k+i].x;
 39           a[j+k].x:=u.x+v.x;
 40           a[j+k].y:=u.y+v.y;
 41           a[j+k+i].x:=u.x-v.x;
 42           a[j+k+i].y:=u.y-v.y;
 43           u:=w;
 44           w.x:=u.x*p.x-u.y*p.y;
 45           w.y:=u.x*p.y+u.y*p.x;
 46         end;
 47         inc(j,s);
 48       end;
 49       i:=i shl 1;
 50     end;
 51     if ty=-1 then
 52     begin
 53       for i:=0 to n-1 do
 54         a[i].x:=a[i].x/extended(n);
 55      end;
 56   end;
 57
 58 begin
 59   readln(n);
 60   m:=n;
 61   n:=1;
 62   while n<=2*(m-1) do
 63   begin
 64     n:=n shl 1;
 65     inc(l);
 66   end;
 67   for i:=0 to n-1 do
 68     r[i]:=(r[i shr 1] shr 1) or ((i and 1) shl (l-1));
 69   for i:=0 to m-1 do
 70     readln(d[i]);
 71   for i:=0 to m-1 do
 72   begin
 73     a[i].x:=d[i];
 74     if i<>0 then b[i].x:=1/i/i;
 75   end;
 76   fft(a,1); fft(b,1);
 77   for i:=0 to n-1 do
 78   begin
 79     c[i].x:=a[i].x*b[i].x-a[i].y*b[i].y;
 80     c[i].y:=a[i].x*b[i].y+a[i].y*b[i].x;
 81   end;
 82   fft(c,-1);
 83   for i:=0 to m-1 do
 84     ans[i]:=c[i].x;
 85
 86   for i:=0 to n-1 do
 87   begin
 88     a[i].x:=0;
 89     a[i].y:=0;
 90     if i<m then a[i].x:=d[m-1-i];
 91   end;
 92   fft(a,1);
 93   for i:=0 to n-1 do
 94   begin
 95     c[i].x:=a[i].x*b[i].x-a[i].y*b[i].y;
 96     c[i].y:=a[i].x*b[i].y+a[i].y*b[i].x;
 97   end;
 98   fft(c,-1);
 99   for i:=0 to m-1 do
100   begin
101     ans[i]:=ans[i]-c[m-i-1].x;
102     writeln(ans[i]:0:3);
103   end;
104 end.

bzoj3160

好题，首先合法方案=不连续间隔相等的回文-连续的回文

连续回文的方案显然可以manacher搞出来

连续的呢？考虑穷举中心点，我们只要快算计算中间点两边的相等间隔字符相等的对数s

则此中心点贡献的方案数即为2^s-1

设中心点为k，则s=sigama([a(i)=a(k-i)]) 这里k代表的是manacher翻倍后的下标，i为原串的下标

考虑字符只有a，b两种，我们只要分别统计即可

当统计a时，把a标为1，b标为0，则[a(i)=a(k-i)]=a(i)*a(k-i) 卷积就出现了，fft搞搞即可

注意inline可以加快operator的速度，但要慎用，有时候会出错

  1 uses math;
  2 const mo=1000000007;
  3 type point=record
  4        x,y:double;
  5      end;
  6      arr=array[0..300010] of point;
  7
  8 var a,b:arr;
  9     p,d,r:array[0..300010] of longint;
 10     c:array[0..300010] of char;
 11     s:ansistring;
 12     x,n,l,len,i,ans:longint;
 13
 14 operator +(a,b:point)c:point;inline;
 15   begin
 16     c.x:=a.x+b.x;
 17     c.y:=a.y+b.y;
 18   end;
 19
 20 operator -(a,b:point)c:point;inline;
 21   begin
 22     c.x:=a.x-b.x;
 23     c.y:=a.y-b.y;
 24   end;
 25
 26 operator *(a,b:point)c:point;inline;
 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 function min(a,b:longint):longint;
 33   begin
 34     if a>b then exit(b) else exit(a);
 35   end;
 36
 37 function manacher:longint;
 38   var w,t,i,right,k:longint;
 39   begin
 40     c[0]:=‘$‘;
 41     t:=0;
 42     for i:=1 to len do
 43     begin
 44       inc(t); c[t]:=‘#‘;
 45       inc(t); c[t]:=s[i];
 46     end;
 47     c[t+1]:=‘#‘;
 48     right:=0; k:=0;
 49     w:=0;
 50     for i:=1 to t do
 51     begin
 52       if right>i then p[i]:=min(p[2*k-i],right-i)
 53       else p[i]:=1;
 54       while c[i+p[i]]=c[i-p[i]] do inc(p[i]);
 55       w:=(w+p[i] div 2) mod mo;
 56       if p[i]+i>right then
 57       begin
 58         right:=p[i]+i;
 59         k:=i;
 60       end;
 61     end;
 62     exit(w);
 63   end;
 64
 65 procedure fft(var a:arr; ty:longint);
 66   var i,j,k,s:longint;
 67       w,p,u,v:point;
 68   begin
 69     for i:=0 to n-1 do
 70       if i<r[i] then
 71       begin
 72         w:=a[r[i]];
 73         a[r[i]]:=a[i];
 74         a[i]:=w;
 75       end;
 76     i:=1;
 77     while i<n do
 78     begin
 79       p.x:=cos(pi/i);
 80       p.y:=ty*sin(pi/i);
 81       s:=i shl 1;
 82       j:=0;
 83       while j<n do
 84       begin
 85         w.x:=1; w.y:=0;
 86         for k:=0 to i-1 do
 87         begin
 88           u:=a[j+k];
 89           v:=w*a[j+k+i];
 90           a[j+k]:=u+v;
 91           a[j+k+i]:=u-v;
 92           w:=w*p;
 93         end;
 94         inc(j,s);
 95       end;
 96       i:=i shl 1;
 97     end;
 98
 99   end;
100 begin
101   readln(s);
102   len:=length(s);
103   d[0]:=1;
104   for i:=1 to len do
105     d[i]:=d[i-1]*2 mod mo;
106   n:=1;
107   l:=0;
108   while n<=(len-1) shl 1 do
109   begin
110     n:=n*2;
111     inc(l);
112   end;
113   for i:=0 to n-1 do
114     r[i]:=(r[i shr 1] shr 1) or ((i and 1) shl (l-1));
115   for i:=1 to len do
116     if s[i]=‘a‘ then a[i-1].x:=1;
117   fft(a,1);
118   for i:=0 to n-1 do
119   begin
120     b[i]:=a[i]*a[i];
121     a[i].x:=0; a[i].y:=0;
122   end;
123   for i:=1 to len do
124     if s[i]=‘b‘ then a[i-1].x:=1;
125   fft(a,1);
126   for i:=0 to n-1 do
127     b[i]:=b[i]+a[i]*a[i];
128   fft(b,-1);
129   for i:=0 to n-1 do
130   begin
131     x:=trunc(round(b[i].x/n));
132     ans:=(ans+d[(x+1) shr 1]-1) mod mo;
133   end;
134   writeln((ans-manacher+mo) mod mo);
135 end.