转自 z55250825 的几篇关于FFT的博文(二)

题目大意:高精度乘法。

fft的实现貌似有很多种,咱先写的是一种递归的fft,应该算是比较快的了吧。参考了 Evil君 的代码,那个运算符重载看的咱P党泪流满面。 (没想到P竟然有运算符重载咩...)

先背模板再理解0.0

以下是待补的对模板的理解

{

其实讲的主要的关键就是如何递归,他记录了一个深度 t,一个左边界s(开区间的),和一个最后FFT的结果的数组a。

他实际上是在递归的过程中就已经计算好了叶子的了,所以复杂度是O(nlogn),咱们来看看咱们如何通过这些递归变量计算出咱们需要的需要计算的(实际上弄出这个来这个算法就基本上可以了)。

主要的代码段就是这里:

{  

 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];

回顾之前的FFT算法,咱们在深度为K的时候,将A分成两半,实际上可以看做是将A按照二进制的第K位是否为1分成两半的。以这个咱们发现..咱们记录的那个S实际上就是当前进行FFT所包含的A的元素的下标二进制的公共部分....然后当前深度的合并的元素实际上就在这里面,且就是在S的二进制基础上在最高位加上一个二进制,即0~ n shr t,即S+[0,n shr t-1]shl (t+1)(这个实际上就是代码中的p了 0v0)

咱们攻克了第二行,可是第一行的那个循环还是有点不对劲的感觉...按道理应该是 0~n shr t-1的,可是这里是0~n shr t shr 1-1,实际上这就是要用蝴蝶操作了,所以只需要枚举一半的量,也就是说,这里实际上要确定两个进行蝴蝶操作的A的下标的二进制关系。

首先确定一点,咱们在整个算法中用到a的地方仅仅在于划分,所以咱们用它们的对应的位置存储y[],然后之后都是直接用 Y0+xY1这样的形式来计算,所以咱们要清楚循环里面存的实际上应该是y。当前所要计算的应该是 w(n>>t,)系列的,咱们将它们存储到对应的区间里面。(其实就是原来的!!!)

对于某一个区间给定的 s,t,咱们首先可以计算出这个里面涉及到的 a[],为[s,s+n>>t],然后分段,一段就是 i+1<<(t+1),其对应的应该就是+1<<t.

}

吐槽:这个模板还需要优化额..为什么感觉和普通的高精度比还是不太行= =

(不过去wikioi的那个10^5的测试貌似跑的还是蛮快的...但是目测常数略大...,因为10^5的数要800ms+跑完)

========================

program fft;

type cp=record x,y:double;end;

     arr=array[0..1 shl 14]of cp;

var a,b,tt,w:array[0..1 shl 14]of cp;

    c:array[0..1000010]of longint;

    n,tot1,tot2:longint;

 

operator *(var a,b:cp)c:cp;

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;

begin c.x:=a.x+b.x;c.y:=a.y+b.y;end;

operator -(var a,b:cp)c:cp;

begin c.x:=a.x-b.x;c.y:=a.y-b.y;end;

 

procedure dft(var a:arr;s,t:longint);

var i,p:longint;

    wt:cp;

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 init;

var ch:char;

    i,k:longint;

    j:cp;

begin

  read(ch);tot1:=0;tot2:=0;

  while (ord(ch)>=ord(‘0‘))and(ord(ch)<=ord(‘9‘)) do

   begin

     a[tot1].x:=ord(ch)-ord(‘0‘);

     read(ch);

     inc(tot1);

   end;

  read(ch);

  while (ord(‘0‘)<=ord(ch)) and(ord(ch)<=ord(‘9‘)) do

   begin

     b[tot2].x:=ord(ch)-ord(‘0‘);

     read(ch);

     inc(tot2);

   end;

  dec(tot1);dec(tot2);

  for i:=0 to tot1 shr 1 do

   begin

     j:=a[i];a[i]:=a[tot1-i];a[tot1-i]:=j;

   end;

  for i:=0 to tot2 shr 1 do

   begin

     j:=b[i];b[i]:=b[tot2-i];b[tot2-i]:=j;

   end;

  if tot1<tot2 then tot1:=tot2;

  n:=1;

  while n>>1<(tot1+1) do n:=n shl 1;

  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 w[i].y:=-w[i].y;

  for i:=0 to n-1 do a[i]:=a[i]*b[i];

  dft(a,0,0);

  fillchar(c,sizeof(c),0);

  for i:=0 to n-1 do

   begin

     c[i]:=c[i]+round(a[i].x/n);

     c[i+1]:=c[i] div 10;

     c[i]:=c[i] mod 10;

   end;

  i:=n;

  while (c[i]=0)and(i>0) do dec(i);

  for k:=i downto 0 do write(c[k]);

end;

 

begin

  init;

 end.

转自 z55250825 的几篇关于FFT的博文(二)

时间: 2024-11-09 02:44:43

转自 z55250825 的几篇关于FFT的博文(二)的相关文章

转自 z55250825 的几篇关于FFT的博文(一)

关于FFT,咱们都会迫不及待地 @  .....(大雾)(貌似被玩坏了...) .....0.0学习FFT前先orz FFT君. 首先先是更详细的链接(手写版题解点赞0v0) FFT的资料 其实众所周知的最详细的算法解释在<算法导论>上...然后咱就是边看着那个边码理解的... 首先来看看多项式乘法和快速FFT的关系,然后咱们再来看能否聊到卷积什么的东西... 其实觉得还是去看算法导论最好. [一.多项式及其表达方式.] 首先什么是多项式额.....实际上就是一个类似这个的东西. A(x)=∑

转自 z55250825 的几篇关于FFT的博文(三)

题目大意:给出n个数qi,定义 Fj为 令 Ei=Fi/qi,求Ei. 其实这道题就是看到有FFT模板才觉得有必要学一下的... 所以实际上就是已经知道题解了... = =. 所以问题就是求这两个多项式相乘的系数. 这里咱卷积不太熟悉,所以咱们来证明一下这个结论显然还是不错的. 首先咱们设 f(x)(k) 表示f(x)的 第k项的系数(就是 x^k-1 那一项) 那么首先了解下卷积,如果f,g是一个序列(这指的是其系数构成的序列),那么卷积S也是一个序列 咱们有: S(k)=∑f(x)(i)*g

React-Native入门指南——第4篇react-native布局实战(二)

React-Native入门指南 github:https://github.com/vczero/react-native-lession React-Native:用JavaScript开发你的原生应用,释放Native的UI体验,体验 Hybird开发效率. 最近一个星期写的文章如下,链接是github page的,其实也可以在系列博客找到相应文章: 第1篇hello react-native 第2篇认识代码结构 第3篇css和布局 第4篇学会react-native布局(一) 第4篇re

Python开发【第十八篇】:MySQL(二)

Python开发[第十八篇]:MySQL(二) 视图 视图是一个虚拟表(非真实存在),其本质是[根据SQL语句获取动态的数据集,并为其命名],用户使用时只需使用[名称]即可获取结果集,并可以将其当作表来使用. SELECT * FROM ( SELECT nid, NAME FROM tb1 WHERE nid > 2 ) AS A WHERE A. NAME > 'alex'; 临时表搜索 1.创建视图 --格式:CREATE VIEW 视图名称 AS SQL语句 CREATE VIEW v

关于Java多线程的线程同步和线程通信的一些小问题(顺便分享几篇质量高的博文)

Java多线程的线程同步和线程通信的一些小问题(顺便分享几篇质量高的博文) 前言:在学习多线程时,遇到了一些问题,这里我将这些问题都分享出来,同时也分享了几篇其他博客主的博客,并且将我个人的理解也分享给大家. 一.对于线程同步和同步锁的理解(注:分享了三篇高质量的博客) 以下我精心的挑选了几篇博文,分别是关于对线程同步的理解和如何选择线程锁以及了解线程锁的作用范围. <一>线程同步锁的选择 1. 这里我推荐下Java代码质量改进之:同步对象的选择这篇博文. 2. 以上推荐的博文是以卖火车票为例

iOS开发——语法篇OC篇&amp;高级语法精讲二

Objective高级语法精讲二 Objective-C是基于C语言加入了面向对象特性和消息转发机制的动态语言,这意味着它不仅需要一个编译器,还需要Runtime系统来动态创建类和对象,进行消息发送和转发.下面通过分析Apple开源的Runtime代码(我使用的版本是objc4-646.tar)来深入理解Objective-C的Runtime机制. Runtime数据结构 在Objective-C中,使用[receiver message]语法并不会马上执行receiver对象的message方

【浅墨Unity3D Shader编程】之七 静谧之秋篇: 表面着色器的写法(二)——自定义光照模式

本系列文章由@浅墨_毛星云 出品,转载请注明出处. 文章链接:http://hpw123.net/plus/view.php?aid=183 作者:毛星云(浅墨)    微博:http://weibo.com/u/1723155442 邮箱: [email protected] QQ交流群:330595914 更多文章尽在:http://www.hpw123.net 本文主要讲解了Unity中SurfaceShader的自定义光照模式的写法. 上篇文章中我们已经说到,表面着色器将分为两次讲解,上

Python开发【第九篇】:HTML (二)

python[第十四篇]HTML基础 时间:2016-08-08 20:57:27      阅读:49      评论:0      收藏:0      [点我收藏+] 什么是HTML? HTML(HyperText MarkUp Language)超文本标记语言,通过使用标记来描述文档结构和表现形式的一种语言,由浏览器进行解析,然后把结果显示在网页上,通俗的讲它就是服务器发送的字符串到浏览器,通过浏览器能解析的规则用HTML来描述, 它是网页构成的基础,你见到的所有网页都离不开HTML,所以

【转】一篇关于迭代器的博文

http://blog.csdn.net/sandy_zc_1/article/details/6529304 CSDN Blog上sandy_zc_1的一篇博文,解答了我关于list迭代器和vector,deque迭代器的困惑,受教!