zoukankan      html  css  js  c++  java
  • 关于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.
    View Code

    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.
    View Code

    bzoj3527

    给一下题面吧,【BZOJ3527】【Zjoi2014】【力】 - z55250825 - z55250825, 令 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.
    View Code

    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.
    View Code
  • 相关阅读:
    UVa 116 单向TSP(多段图最短路)
    POJ 1328 Radar Installation(贪心)
    POJ 1260 Pearls
    POJ 1836 Alignment
    POJ 3267 The Cow Lexicon
    UVa 1620 懒惰的苏珊(逆序数)
    POJ 1018 Communication System(DP)
    UVa 1347 旅行
    UVa 437 巴比伦塔
    UVa 1025 城市里的间谍
  • 原文地址:https://www.cnblogs.com/phile/p/4665374.html
Copyright © 2011-2022 走看看