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
  • 相关阅读:
    012_DRC检查与处理
    深度系统20.3中亿图图示任务栏名称显示乱码
    deepin20.3+nvidia460.27+cuda11.2+cudnn8.1.1+anconda3.2021.11+paddle2.1.2
    C++中使用DOM写XML文档
    理解lvalue和rvalue
    C++/CLI与C#常用语法对比
    VC++ MSXML创建XML文件以及对XML文档解析
    Stack overflow 编译能通过,运行时出现Stack overflow
    于typedef的用法总结
    VC2008操作Excel2007总结
  • 原文地址:https://www.cnblogs.com/phile/p/4665374.html
Copyright © 2011-2022 走看看