[关闭]
@zzzc18 2018-04-17T11:06:23.000000Z 字数 37375 阅读 3413

快速傅里叶变换学习笔记

FFT bzoj

本篇博客中题目的代码均位于文末，文中不再出现

算法原理

Miskcoo的讲解

对于快速傅里叶变换的原理，可以参考上面的博客以及《算法导论》

大体说明

FFT用于加速多项式的乘法：


即：

上面两个式子都是卷积的形式，也就是FFT可以优化的类型

使用暴力做法，时间复杂度为 ，而使用快速傅里叶变换可以做到，于是许多问题从不可做变成了可做，FFT对生成函数的优化是个很好的例子，可以快速解决组合问题。

题目

模板

UOJ#34：多项式乘法

CODE

或者bzoj2179

FFT优化高精度乘法，一个数就可以看作，其中：

BZOJ2194: 快速傅立叶之二

Description

请计算其中 ，并且有 。 , 中的元素均为小于等于100的非负整数。 , 最高次项均为

Solution

没错，就是一个多项式的乘法

我们将上面那个式子变一下形
首先将 的系数全部 reverse
eg:


此时

然后改变求和指标，上式变为

另 ，则 ，代入上式可得

与之前类似地：

因此我们只需要将 反转，得到的卷积结果是 的反转

BZOJ3527: [Zjoi2014]力

Description

给出 个数 ，给出 的定义如下：

令 ，求 .

Solution

很显然，

另 ，，那么：

而另一部分与BZOJ2194类似：

两次答案加起来就是最终答案

另：这题卡精度，计算 时要 ，不能

BZOJ3771: Triple

Description

我们讲一个悲伤的故事。
从前有一个贫穷的樵夫在河边砍柴。
这时候河里出现了一个水神，夺过了他的斧头，说：
“这把斧头，是不是你的？”
樵夫一看：“是啊是啊！”
水神把斧头扔在一边，又拿起一个东西问：
“这把斧头，是不是你的？”
樵夫看不清楚，但又怕真的是自己的斧头，只好又答：“是啊是啊！”
水神又把手上的东西扔在一边，拿起第三个东西问：
“这把斧头，是不是你的？”
樵夫还是看不清楚，但是他觉得再这样下去他就没法砍柴了。
于是他又一次答：“是啊是啊！真的是！”
水神看着他，哈哈大笑道：
“你看看你现在的样子，真是丑陋！”
之后就消失了。
樵夫觉得很坑爹，他今天不仅没有砍到柴，还丢了一把斧头给那个水神。
于是他准备回家换一把斧头。
回家之后他才发现真正坑爹的事情才刚开始。
水神拿着的的确是他的斧头。
但是不一定是他拿出去的那把，还有可能是水神不知道怎么偷偷从他家里拿走的。
换句话说，水神可能拿走了他的一把，两把或者三把斧头。
樵夫觉得今天真是倒霉透了，但不管怎么样日子还得过。
他想统计他的损失。
樵夫的每一把斧头都有一个价值，不同斧头的价值不同。总损失就是丢掉的斧头价值和。
他想对于每个可能的总损失，计算有几种可能的方案。
注意：如果水神拿走了两把斧头 和 ， 和 视为一种方案。拿走三把斧头时， 视为一种方案。

Description2

嗯，一句话，有 个物品各有价值，要求从中取出 到 个物品，不同排列视作同一种，求可能的取出物品的价值和有多少种

所有数据满足：

Solution

我们直接构造出这些斧头的生成函数

举个例子，有斧子们

则其对应的生成函数为

在这一题里，我们记取一把斧子的生成函数为 ，取两把相同的斧子的生成函数为 ， 三把为 那么取一个斧头的方案就是 ；
对于两个斧头，方案数变为了
对于三把斧头，方案数就是
上边的计算，只要理解生成函数再使用容斥原理就可以了

注：对于函数的操作，可以直接转化为点值表达后一起搞事，不用每次 、，可以参见代码

BZOJ3509: [CodeChef] COUNTARI

Description

给定一个长度为 的数组 ，求有多少对 （满足 。

Solution

首先，通过移项，不难发现，我们要求的是对于每一个 满足 的 对数 （
所以我们直接从头到尾扫一遍，分别维护前缀及后缀里各个值出现的次数，看对于当前 有那些组合满足条件。时间复杂度

此时我们再来看这个式子，这个也是可以使用生成函数的。
先把前缀及后缀的生成函数搞出来，分别记为 。
令 ，那么，函数 的第 项的系数就是对于当前 满足条件的对数，卷积可以用FFT优化到 。
可惜的是这玩意最终的时间复杂度是 ，硬上并无卵用

用了FFT还是这么慢的主要原因是用了太多次，可不可以用少一些呢？
现在我们将待处理的数列分为 个块，对于块内，暴力解决，而对于块外则使用FFT。
具体来说就是在一个块 中，可以通过暴力统计计算出满足 以及 的方案数；
而对于，我们构造 的生成函数以及 的生成函数，卷积一下就是方案数了；
最终两个答案加一起就是最终答案

另：这题要注意常数

BZOJ3456: 城市规划

Description

刚刚解决完电力网络的问题, 阿狸又被领导的任务给难住了.
刚才说过, 阿狸的国家有 个城市,现在国家需要在某些城市对之间建立一些贸易路线,使得整个国家的任意两个城市都直接或间接的连通.为了省钱,每两个城市之间最多只能有一条直接的贸易路径.对于两个建立路线的方案, 如果存在一个城市对, 在两个方案中是否建立路线不一样, 那么这两个方案就是不同的, 否则就是相同的. 现在你需要求出一共有多少不同的方案.好了, 这就是困扰阿狸的问题. 换句话说, 你需要求出 个点的简单(无重边无自环)无向连通图数目.由于这个数字可能非常大,你只需要输出方案数 即可.

Solution

首先对于这个模数，已经暴露了题目想让我们搞什么——快速数论变换

下文直接用 表示

我们记最终的答案为一个函数 ，那么 是最终答案。
如果不考虑这个图的连通性，那么答案就显然是 ，即所有可能的边选或者不选。设
进一步考虑，如何找出 与 的等量关系？
我们枚举 号点所在的连通块，我们可以得到以下的等量关系：


同时除以 得

是不是看到了卷积？美滋滋

令：（这样写好像不规范）


其中为了求逆元


只要求一个 逆元一卷积就好了（求逆元的方法回头再说）

BZOJ4503: 两个串

Description

兔子们在玩两个串的游戏。给定两个字符串 和 ，兔子们想知道 在 中出现了几次，
分别在哪些位置出现。注意T中可能有 字符，这个字符可以匹配任何字符。
长度不超过 ， 长度不会超过 。 中只包含小写字母，T中只包含小写字母和

Solution

如果两个长度为 的串 相等时，我们可以写成

也可以写成

现在将 翻转，得到 ，上式就可以写成

我们记 中所对应的各字符减去 的值记为 ，其中当字符为 时，将这个值赋为 ，就可以有一下关系( 的系数颠倒后为 )

感觉不用解释了， 的某一位是 的时候就意味着有一个满足题意的匹配

BZOJ3160: 万径人踪灭

Description

T1

T2

T3

T4

Solution

并没有看完整题意。。。

首先，如果回文串要求连续的，可以通过 算法 内很轻松解决
这个Manacher讲的不错
那此时就可以计算所有的，不管连续不连续的回文串，再减去连续的就行了

对于计算回文串个数，可以通过一下式子

解释一下，由于只有 两种字符，我们把数列中的 赋为 , 赋为 ，做一次计算；
再把数列中的 赋为 , 赋为 ，做一次计算；两次结果相加。
这样可以得到以 为中心的、最长的、不管连续不连续的回文串。
那么其数量为 。
最后把所有答案加起来减去连续的回文串数量就行

BZOJ3992: [SDOI2015]序列统计

Description

小C有一个集合 ，里面的元素都是小于 的非负整数。他用程序编写了一个数列生成器，可以生成一个长度为N的数列，数列中的每个数都属于集合 。
小C用这个生成器生成了许多这样的数列。但是小C有一个问题需要你的帮助：给定整数 ，求所有可以生成出的，且满足数列中所有数的乘积的值等于 的不同的数列的有多少个。小C认为，两个数列和不同，当且仅当至少存在一个整数 ，满足 。另外，小C认为这个问题的答案可能很大，因此他只需要你帮助他求出答案的值就可以了。

输入：一行，四个整数，N、M、x、|S|，其中|S|为集合S中元素个数。第二行，|S|个整数，表示集合S中的所有元素。

对于10%的数据，1<=N<=1000；

对于30%的数据，3<=M<=100；

对于60%的数据，3<=M<=800；

对于全部的数据，1<=N<=109，3<=M<=8000，M为质数，1<=x<=M-1，输入数据保证集合S中元素不重复

Solution

模数又暴露了算法，但并不影响这是个神题

首先 是个质数，我们可以找到其原根
原根是什么？

假设一个数 对于 来说是原根，那么 的结果两两不同,且有 ,那么 可以称为是 的一个原根
简单来说， （ 为素数）

关于其计算，参见这里

简而言之，每次枚举原根，暴力看是否满足条件

  1. //qpow(x,k,MOD)为快速幂 : x^k %MOD
  2. bool judge(int val){
  3.     if(m==2)
  4.         return 1;
  5.     if(val==1)
  6.         return 0;
  7.     for(int i=2;i*i<=m;i++)
  8.         if((m-1)%i==0)
  9.             if(qpow(val,m/i,m)==1)
  10.                 return 0;
  11.     return 1;
  12. }
  14. int FindG(int val){
  15.     int re=1;
  16.     while(!judge(re)) re++;
  17.     return re;
  18. }

另 的原根是
这时，每个数 就可以写成
所以

就变成了

这一步将乘法运算变为了加法

现在将 构造生成函数记为
显然，因为每个数可以取多次


的第 位就是最终答案

时间复杂度

注：
1.多项式快速幂的时候不能把次数高于 的项直接清零，而是应该将第 项加到 项上
2.对于 特别计算

BZOJ3557: [Ctsc2014]随机数

Description

露露、花花和萱萱最近对计算机中的随机数产生了兴趣。为大家所熟知的是，有计算机生成的随机数序列并非真正的随机数，而是由一定法则生成的伪随机数。 某一天，露露了解了一种生成随机数的方法，成为Mersenne twister。给定初始参数 ，和初值 ，它通过下列递推式构造伪随机数列:
公式

其中XOR是二进制异或运算（C/C++中的^运算）。而参数x的选取若使得该数列在长度趋于无穷时，近似等概率地在 中取值，
就称x为好的。例如，在 时 就显然不是好的。
在露露向伙伴们介绍了Mersenne twister之后，花花想用这一些经典的随机性测试来检验它的随机性强度。为此，花花使用计算机计算
了一些 。
但细心的萱萱注意到，花花在某次使用二进制输入 时，在末尾多输入了 个 。她正想告诉花花这个疏忽，然而花花已经计算并记录了错误的 而没有记录 的值。虽然这其实不是什么致命的问题，但是在萱萱告诉花花她这个疏漏时，作为完美主义者的花花还是恳求萱萱帮她修正 的值。萱萱便把这个任务交给了她的AI——你。

Input:
第一行包含一个正整数 ，
第二行为二进制表示的 （共 个数，从低位到高位排列）
第三行为二进制表示的 （排列方式同 ），
第四行包含一个整数 。
接下来分为两种可能的情况：
1. （萱萱记下了花花的输入）：则第五行包含一个整数，表示萱萱记下来的正确的k值。
2. （萱萱未能记下花花的输入）：则第五行为 ，第六行输入花花计算出错误的二进制表示的 。

Output:
仅一行，为 位的 串，表示你求得的正确 （同样要求从低位到高位）。

Solution

我们可以直接将 看做两个 （系数mod2）下的多项式
在此意义下( 为异或)


所以上面两个操作看成

以及

这样两种操作都可以写成

仔细分析上式是对的
1.操作一由于最高次项不到 取模是没有影响的
2.操作二就是这么个东西

于是 ，可以直接计算出

需要多项式快速幂以及多项式取模（和除法在一起）

对于
因为题目保证 是“好的”，每次操作都是确定的，可以发现一个性质， 在 里恰好每个值出现过一遍，否则 不可能是好的，最后不是等概率。
将上面这个性质写成式子就是

我们再将问题表示出来

再利用上面的性质整理

所以

所以

所以只要求出 在 下的逆元，与 乘起来计算快速幂，再乘上 就行了

求出 在 下的逆元需要拓展欧几里得，这与数论里的方法类似

吐槽一下：
1.这题真是个板子大全（没有开根）
2.真不可调试//用了一整个上午，还是有标程的情况下，代码里的DEBUG（）可以看出我的绝望
3.原题一个点时限20sec,搞到测试数据一测最大点15sec，然而BZOJ丧心病狂总时间160sec(20个点)，我现在已经没心情去卡常了，就这了（这意味着我的代码会TLE，但是肯定对）
4.%%%WJMZBMR，%%%miskcoo

Code

BZOJ3557

  1. /***************************************************************
  2.     Problem: 3557
  3.     User: zzzc18
  4.     Language: C++
  5.     Result: Accepted
  6.     Time:157680 ms
  7.     Memory:167432 kb
  8. ****************************************************************/
  10. #pragma GCC target("mmx,avx,sse4,sse3,sse2,tune=native")
  11. #include <algorithm>
  12. #include   <cstring>
  13. #include    <cctype>
  14. #include    <cstdio>
  15. #include     <cmath>
  16. using namespace std;
  18. namespace zzzc18{
  19.     const int L = (1<<23)+1;
  20.     char ibuf[L],obuf[L];
  21.     char *iS=ibuf,*iT=ibuf,*oS=obuf,*oT=obuf+L-1;
  22. #define Rdc() (iS==iT?(iT=(iS=ibuf)+fread(ibuf,1,L,stdin),iS==iT?EOF:*iS++):*iS++)
  24.     void flush(){
  25.         fwrite(obuf,1,oS-obuf,stdout);
  26.         oS=obuf;
  27.     }
  29.     void Ptc(char ch){
  30.         if(oS==oT)flush();
  31.         *oS++=ch;
  32.     }
  34.     template<typename type>
  35.         void Rd(type &x){
  36.             x=0;
  37.             char ch=Rdc();
  38.             while(ch<'0' || ch>'9'){ch=Rdc();}
  39.             while(ch>='0' && ch<='9'){x=x*10+(ch^'0');ch=Rdc();}
  40.         }
  42.     template<typename type>
  43.         void Pt(type x,char ch=0){
  44.             static char tmp[20];
  45.             static int top;
  46.             if(!x)Ptc('0');
  47.             while(x){tmp[++top]=x%10+'0';x/=10;}
  48.             while(top){Ptc(tmp[top--]);}
  49.             if(ch)
  50.                 Ptc(ch);
  51.         }
  52. }
  53. using zzzc18::Rd;
  54. using zzzc18::Ptc;
  56. typedef long long LL;
  57. const int MAXN = 2100000+9;
  58. const int MOD = 479<<21|1;
  59. const int G = 3;
  60. int X[MAXN],M0[MAXN],tp[MAXN],A1[MAXN],B1[MAXN],A3[MAXN];
  61. int size,bit_length;
  62. int loc[MAXN];
  63. int m;
  64. int tmp[5][MAXN];
  65. int type;
  66. int tot_len;
  67. int A2[MAXN],B2[MAXN];
  68. int W[MAXN],inv_W[MAXN];
  69. int NTT_TOT;
  70. struct PAIR{
  71.     int first,second;
  72.     PAIR(int _x=0,int _y=0):first(_x),second(_y){}
  73. };
  75. void DEBUG(int *A,int len=32){
  76.     for(int i=0;i<len;i++)
  77.         printf("%d ",A[i]);
  78.     puts("");
  79. }
  81. int qpow(int x,int k){
  82.     int re=1;
  83.     int mul=x;
  84.     for(register int i=1;i<=k;i<<=1){
  85.         if(i&k)
  86.             re=(LL)re*mul%MOD;
  87.         mul=(LL)mul*mul%MOD;
  88.     }
  89.     return re;
  90. }
  92. inline void GF2(int &x){
  93.     if(x>=3000000) x=(MOD-x)&1;//x由于减法而+过MOD，所以这样 
  94.     else x&=1;//mod 2
  95. }
  97. int inv(int x){
  98.     return qpow(x,MOD-2);
  99. }
  101. inline void init(int x){
  102.     for(size=1,bit_length=-1;size<x;size<<=1,bit_length++);
  103.     loc[0]=0;
  104.     for(int i=1;i<size;i++){
  105.         loc[i]=loc[i>>1]>>1|((i&1)<<bit_length);
  106.     }
  107. }
  109. void init_eps(int o){
  110.     NTT_TOT=o;
  111.     W[0]=inv_W[0]=1;
  112.     int base=qpow(G,(MOD-1)/o);
  113.     int invbase=inv(base);
  114.     for(int i=1;i<o;i++){
  115.         W[i]=1LL*W[i-1]*base%MOD;
  116.         inv_W[i]=1LL*inv_W[i-1]*invbase%MOD;
  117.     }
  118. }
  120. inline int dec(int x,int y){
  121.     x-=y;
  122.     return x<0?x+MOD:x;
  123. }
  125. inline int inc(int x,int y){
  126.     x+=y;
  127.     return x<MOD?x:x-MOD;
  128. }
  130. void DFT(int *A,int *re){
  131.     for(register int i=0;i<size;i++)
  132.         re[i]=A[loc[i]];
  133.     for(register int k=2;k<=size;k<<=1){
  134.         int len=k>>1;
  135.         int r=NTT_TOT/k;
  136.         for(register int i=0;i<size;i+=k){
  137.             int *x=re+i;int *y=re+i+len;
  138.             for(register int j=0,s=0;j<len;j++,s+=r){
  139.                 int u=x[j],v=(LL)W[s]*y[j]%MOD;
  140.                 x[j]=inc(u,v);
  141.                 y[j]=dec(u,v);
  142.             }
  143.         }
  144.     }
  145. }
  147. void IDFT(int *A,int *re){
  148.     for(register int i=0;i<size;i++)
  149.         re[i]=A[loc[i]];
  150.     for(register int k=2;k<=size;k<<=1){
  151.         int len=k>>1;
  152.         int r=NTT_TOT/k;
  153.         for(register int i=0;i<size;i+=k){
  154.             int *x=re+i;int *y=re+i+len;
  155.             for(register int j=0,s=0;j<len;j++,s+=r){
  156.                 int u=x[j],v=(LL)inv_W[s]*y[j]%MOD;
  157.                 x[j]=inc(u,v);
  158.                 y[j]=dec(u,v);
  159.             }
  160.         }
  161.     }
  162.     int tmp=inv(size);
  163.     for(int i=0;i<size;i++)
  164.         re[i]=(LL)re[i]*tmp%MOD;
  165. }
  167. int GetInv(int degree,int *A,int *B){
  168.     if(degree==1){
  169.         B[0]=inv(A[0]);
  170.         return degree;
  171.     }
  172.     int val=GetInv((degree+1)>>1,A,B);
  173.     init(degree<<1);
  175.     copy(A,A+degree,tp);
  176.     //memcpy(tp,A,(degree-1)<<2);
  177.     fill(tp+degree,tp+size,0);
  178.     DFT(tp,A2);
  180.     fill(B+val,B+size,0);
  181.     DFT(B,B2);
  183.     for(register int i=0;i<size;i++)
  184.         A2[i]=(LL)A2[i]*B2[i]%MOD;
  185.     IDFT(A2,A3);
  186.     for(register int i=0;i<size;i++)
  187.         GF2(A3[i]);
  188.     DFT(A3,A2);
  189.     for(register int i=0;i<size;i++)
  190.         B2[i]=(LL)B2[i]*dec(2,A2[i])%MOD;
  191.     IDFT(B2,B);
  192.     for(register int i=0;i<degree;i++)
  193.         GF2(B[i]);
  194.     fill(B+degree,B+size,0);
  195.     return degree;
  196. }
  198. // A(x)=D(x)B(x)+R(x)
  199. void DivMod(int len_a,int len_b,int *A,int *B,int *D,int *R){
  200.     static int tp[MAXN];
  201.     static int Z[MAXN];
  202.     static bool solved=0;
  203.     int t=len_a-len_b+1;
  204.     init(t<<1);
  205.     if(!solved || type){
  206.         if(!type)solved=1;////对于type==0，计算时模数均为X，只算一遍
  207.         fill(B1,B1+size,0);
  208.         reverse_copy(B,B+len_b,B1);
  209.         GetInv(t,B1,A1);
  210.         fill(A1+t,A1+size,0);
  211.         for(register int i=0;i<t;i++) GF2(A1[i]);
  212.         DFT(A1,Z);
  213.     }
  215.     reverse_copy(A,A+len_a,tp);
  216.     fill(tp+t,tp+size,0);
  217.     DFT(tp,A1);
  218.     for(register int i=0;i<size;i++)
  219.         tp[i]=(LL)A1[i]*Z[i]%MOD;
  220.     IDFT(tp,A1);
  221.     reverse(A1,A1+t);
  222.     for(register int i=0;i<t;i++) GF2(A1[i]);
  223.     if(D)
  224.         copy(A1,A1+t,D);
  225.     //A1 carrys D
  226.     if(!R)return;
  227.     init(len_a);
  228.     if(t<size)fill(A1+t,A1+size,0);
  229.     DFT(A1,tp);
  230.     DFT(B,B1);
  231.     for(register int i=0;i<size;i++)
  232.         tp[i]=(LL)tp[i]*B1[i]%MOD;
  233.     IDFT(tp,A1);
  234.     for(register int i=0;i<len_b;i++){
  235.         R[i]=dec(A[i],A1[i]);//(A[i]-A1[i]+MOD)%MOD;
  236.         GF2(R[i]);
  237.     }
  238.     fill(R+len_b,R+size,0);
  239.     //DEBUG(A,32);
  240. }
  242. void mul(int len,int *A,int *B,int *C){
  243.     init(len);
  244.     if(A!=B){//用来优化常数
  245.         int pos1,pos2;
  246.         for(pos1=size-1;pos1 && !A[pos1];pos1--);//用来优化常数
  247.         for(pos2=size-1;pos2 && !B[pos2];pos2--);//用来优化常数
  248.         init(pos1+pos2+1);
  249.         DFT(A,A1);DFT(B,B1);
  250.         for(register int i=0;i<size;i++)
  251.             A1[i]=(LL)A1[i]*B1[i]%MOD;
  252.         IDFT(A1,C);
  253.     }
  254.     else{
  255.         int pos1;
  256.         for(pos1=size-1;pos1 && !A[pos1];pos1--);//用来优化常数
  257.         init(pos1+pos1+1);
  258.         //printf("%d %d\n",pos1,size);
  259.         DFT(A,A1);
  260.         for(register int i=0;i<size;i++)
  261.             A1[i]=(LL)A1[i]*A1[i]%MOD;
  262.         IDFT(A1,C);
  263.     }
  264.     //printf("NO1:");DEBUG(C,32);
  265.     for(int i=0;i<size;i++) GF2(C[i]);
  266.     if(size<len)
  267.         fill(C+size,C+len,0);
  268.     bool flag=0;
  269.     for(register int i=m;i<size;i++){
  270.         if(C[i]){
  271.             flag=1;
  272.             break;
  273.         }
  274.     }
  275.     if(flag)//取模对其有影响
  276.         DivMod(len,tot_len,C,X,0,C);//取模
  277.     //printf("NO2:");DEBUG(C,32);
  278. }
  280. PAIR Extgcd(int len_a,int len_b,int *A,int *B,int *XX,int *Y,int *R){
  281.     while(len_a && !A[len_a-1]) len_a--;
  282.     while(len_b && !B[len_b-1]) len_b--;
  283.     if(!len_b){
  284.         XX[0]=1,Y[0]=0;
  285.         return PAIR(1,1);
  286.     }
  287.     if(len_a<len_b){
  288.         PAIR re=Extgcd(len_b,len_a,B,A,Y,XX,R);
  289.         swap(re.first,re.second);
  290.         return re;
  291.     }
  292.     int t=len_a-len_b+1;
  293.     init(t+len_b);
  294.     int *D=new int[size<<1];
  296.     DivMod(len_a,len_b,A,B,D,R);
  297.     //DEBUG(D);
  298.     PAIR len=Extgcd(len_b,len_b,B,R,XX,Y,A);
  299.     /*printf("%d %d\n",len.first,len.second);
  300.       printf("X:");DEBUG(XX,32);
  301.       printf("Y:");DEBUG(Y,32);*/
  302.     init(t+len.second);
  304.     PAIR re;
  305.     copy(XX,XX+len.first,A);
  306.     fill(A+len.first,A+size,0);
  307.     copy(Y,Y+len.second,XX);
  308.     re.first=len.second;
  309.     fill(D+t,D+size,0);fill(Y+len.second,Y+size,0);
  310.     DFT(D,A1);DFT(Y,B1);
  311.     for(int i=0;i<size;i++)
  312.         tp[i]=(LL)A1[i]*B1[i]%MOD;
  313.     delete[] D;
  314.     IDFT(tp,Y);
  315.     for(register int i=0;i<size;i++){
  316.         Y[i]=dec(A[i],Y[i]);
  317.         GF2(Y[i]);
  318.     }
  320.     int l;
  321.     for(l=size;l && !Y[l-1];l--);
  322.     re.second=l;
  323.     return re;
  324. }
  326. void solve0(){
  327.     int k;
  328.     Rd(k);
  329.     tot_len=m+1;
  330.     init(tot_len<<1);
  331.     init_eps(size);
  332.     int lena=0,lenb=1;
  333.     bool flag=0;
  334.     int mdiv2=m>>1;
  335.     int *val=tmp[0];
  336.     int *ans=tmp[1];
  337.     for(register int i=1;i<=k;i<<=1){
  338.         if(flag){
  339.             if(i&k)
  340.                 mul(tot_len<<1,ans,val,ans);
  341.             mul(tot_len<<1,val,val,val);
  342.         }
  343.         else{
  344.             if(i&k)lena+=lenb;
  345.             lenb<<=1;
  346.             if(lena>mdiv2 || lenb>mdiv2){
  347.                 flag=1;
  348.                 val[lenb]=1;
  349.                 ans[lena]=1;
  350.             }
  351.         }
  352.         //printf("val:");DEBUG(val,32);
  353.         //rintf("ans:");DEBUG(ans,32);
  354.     }
  355.     if(!flag)ans[lena]=1;
  356.     mul(tot_len<<1,ans,M0,M0);
  357.     for(int i=0;i<m;i++)
  358.         Ptc('0'+M0[i]);
  359.     Ptc('\n');
  360. }
  362. void solve1(){
  363.     int l;
  364.     Rd(l);
  365.     int *Mk=tmp[0];int *P=tmp[1];
  366.     int *Xval=tmp[2],*Y=tmp[3];
  367.     tot_len=m+1;
  368.     init(max(2<<l,tot_len<<1));
  369.     init_eps(size<<1);
  370.     copy(X,X+size,P);
  371.     copy(M0,M0+size,Mk);
  373.     PAIR len=Extgcd(m,tot_len,Mk,P,Xval,Y,tmp[4]);
  374.     //printf("X:");DEBUG(Xval,32);
  375.     //printf("Y:");DEBUG(Y,32);
  377.     init(max(2<<l,tot_len<<1));
  378.     for(register int i=0;i<m;i++)
  379.         Rd(Mk[i]);
  380.     fill(Mk+m,Mk+size,0);
  381.     //memset(Mk+m,0,(size-m+1)<<2);
  382.     mul(len.first+m,Mk,Xval,Y);
  383.     type=0;
  384.     for(int i=0;i<m-l;i++)
  385.         mul(tot_len<<1,Y,Y,Y);
  386.     mul(tot_len<<1,Y,M0,Y);
  387.     for(int i=0;i<m;i++)
  388.         Ptc('0'+Y[i]);
  389.     Ptc('\n');
  390. }
  392. int main(){
  393.     Rd(m);
  394.     for(register int i=0;i<m;i++)
  395.         Rd(X[i]);
  396.     X[m]=1;
  397.     for(register int i=0;i<m;i++)
  398.         Rd(M0[i]);
  399.     Rd(type);
  400.     if(type)solve1();
  401.     else solve0();
  402.     zzzc18::flush();
  403.     return 0;
  404. }

BZOJ3992

  1. /**************************************************************
  2.     Problem: 3992
  3.     User: zzzc18
  4.     Language: C++
  5.     Result: Accepted
  6.     Time:3180 ms
  7.     Memory:7852 kb
  8. ****************************************************************/
  10. #include<cstdio>
  11. #include<cstring>
  12. #include<algorithm>
  13. using namespace std;
  15. const int MAXM = 300000;
  16. int n,m,X,S;
  17. const int G=3;
  18. const int MOD = (479<<21)|1;
  19. int size,bit_length;
  20. int loc[MAXM];
  21. int num[MAXM];
  22. int fac[MAXM];
  24. int qpow(int x,int k,int mod){
  25.     int re=1;
  26.     int mul=x;
  27.     for(int i=1;i<=k;i<<=1){
  28.         if(i&k)
  29.             re=1LL*re*mul%mod;
  30.         mul=1LL*mul*mul%mod;
  31.     }
  32.     return re;
  33. }
  35. int inv(int x){
  36.     return qpow(x,MOD-2,MOD);
  37. }
  39. void init(int val){
  40.     for(size=1,bit_length=-1;size<val;size<<=1,bit_length++);
  41.     int now=0;
  42.     for(int i=0;i<size;i++){
  43.         loc[i]=now;
  44.         for(int j=1<<bit_length;(now^=j)<j;j>>=1);
  45.     }
  46. }
  48. void DFT(int *A,int *re){
  49.     for(int i=0;i<size;i++)
  50.         re[i]=A[loc[i]];
  51.     for(int k=2;k<=size;k<<=1){
  52.         int len=k>>1;
  53.         int Wn=qpow(G,(MOD-1)/k,MOD);
  54.         for(int i=0;i<size;i+=k){
  55.             int W=1;
  56.             for(int j=0;j<len;j++){
  57.                 int u=re[i+j],v=1LL*W*re[i+j+len]%MOD;
  58.                 re[i+j]=(u+v)%MOD;
  59.                 re[i+j+len]=(u-v+MOD)%MOD;
  60.                 W=1LL*W*Wn%MOD;
  61.             }
  62.         }
  63.     }
  64. }
  66. void IDFT(int *A,int *re){
  67.     for(int i=0;i<size;i++)
  68.         re[i]=A[loc[i]];
  69.     for(int k=2;k<=size;k<<=1){
  70.         int len=k>>1;
  71.         int Wn=inv(qpow(G,(MOD-1)/k,MOD));
  72.         for(int i=0;i<size;i+=k){
  73.             int W=1;
  74.             for(int j=0;j<len;j++){
  75.                 int u=re[i+j],v=1LL*W*re[i+j+len]%MOD;
  76.                 re[i+j]=(u+v)%MOD;
  77.                 re[i+j+len]=(u-v+MOD)%MOD;
  78.                 W=1LL*W*Wn%MOD;
  79.             }
  80.         }
  81.     }
  82.     int tmp=inv(size);
  83.     for(int i=0;i<size;i++)
  84.         re[i]=1LL*re[i]*tmp%MOD;
  85. }
  87. bool judge(int val){
  88.     if(m==2)
  89.         return 1;
  90.     if(val==1)
  91.         return 0;
  92.     for(int i=2;i*i<=m;i++)
  93.         if((m-1)%i==0)
  94.             if(qpow(val,m/i,m)==1)
  95.                 return 0;
  96.     return 1;
  97. }
  99. int FindG(int val){
  100.     static int tmp[MAXM];
  101.     int re=1;
  102.     while(!judge(re)) re++;
  103.     return re;
  104. }
  106. void PRE(){
  107.     int val=FindG(m);
  108.     int now=val;
  109.     for(int i=1;i<m-1;i++){
  110.         fac[now]=i;
  111.         now=1LL*now*val%m;
  112.     }
  113. }
  115. void transform(int *A,int *B,int *ans){
  116.     static int A1[MAXM],B1[MAXM];
  117.     if(A==B){
  118.         DFT(A,A1);
  119.         for(int i=0;i<size;i++)
  120.             A1[i]=1LL*A1[i]*A1[i]%MOD;
  121.         IDFT(A1,ans);
  122.     }
  123.     else{
  124.         DFT(A,A1);DFT(B,B1);
  125.         for(int i=0;i<size;i++)
  126.             A1[i]=1LL*A1[i]*B1[i]%MOD;
  127.         IDFT(A1,ans);
  128.     }
  129.     for(int i=m-1;i<size;i++){
  130.         ans[i%(m-1)]+=ans[i];
  131.         ans[i%(m-1)]%=MOD;
  132.         ans[i]=0;
  133.     }
  134. }
  136. void QPOW(int *A,int k){
  137.     static int tmp[MAXM];
  138.     copy(A,A+m+1,tmp);
  139.     fill(A,A+m+1,0);
  140.     A[0]=1;
  141.     for(int i=1;i<=k;i<<=1){
  142.         if(i&k){
  143.             transform(A,tmp,A);
  144.         }
  145.         transform(tmp,tmp,tmp);
  146.     }
  147. }
  149. int main(){
  150.     scanf("%d%d%d%d",&n,&m,&X,&S);
  151.     if(X==0){
  152.         bool flag=0;
  153.         int x;
  154.         for(int i=1;i<=S;i++){
  155.             scanf("%d",&x);
  156.             if(!x){
  157.                 flag=1;
  158.                 break;
  159.             }
  160.         }
  161.         if(flag)
  162.             printf("%d\n",(qpow(S,n,MOD)-qpow(S-1,n,MOD)+MOD)%MOD);
  163.         else
  164.             printf("0\n");
  165.         return 0;
  166.     }
  167.     PRE();
  168.     init(m<<1|1);
  169.     for(int i=1;i<=S;i++){
  170.         int x;
  171.         scanf("%d",&x);
  172.         if(x)
  173.             num[fac[x]]++;
  174.     }
  175.     QPOW(num,n);
  176.     printf("%d\n",num[fac[X]]);
  177.     return 0;
  178. }

BZOJ4503

  1. /**************************************************************
  2.     Problem: 4503
  3.     User: zzzc18
  4.     Language: C++
  5.     Result: Accepted
  6.     Time:2736 ms
  7.     Memory:27188 kb
  8. ****************************************************************/
  10. #include<cmath>
  11. #include<cstdio>
  12. #include<cstring>
  13. #include<algorithm>
  14. using namespace std;
  16. const int MAXN = 300000+9;
  17. const double PI = acos(-1.0);
  18. struct C{
  19.     double x,y;
  20.     C(double _x=0,double _y=0):x(_x),y(_y){}
  21. };
  22. C operator + (const C &A,const C &B){
  23.     return C(A.x+B.x,A.y+B.y);
  24. }
  25. C operator - (const C &A,const C &B){
  26.     return C(A.x-B.x,A.y-B.y);
  27. }
  28. C operator * (const C &A,const C &B){
  29.     return C(A.x*B.x-A.y*B.y,A.x*B.y+A.y*B.x);
  30. }
  31. C operator * (const double &v,const C &B){
  32.     return C(v*B.x,v*B.y);
  33. }
  35. int size,bit_length;
  36. int loc[MAXN];
  37. char s1[MAXN],s2[MAXN];
  38. int num1[MAXN],num2[MAXN];
  39. int num3[MAXN],num4[MAXN];
  40. int len1,len2;
  41. int ans[MAXN];
  42. C A1[MAXN],B1[MAXN];
  43. C A2[MAXN],B2[MAXN];
  44. int B3;
  46. void init(int x){
  47.     for(size=1,bit_length=-1;size<x;size<<=1,bit_length++);
  48.     int now=0;
  49.     for(int i=0;i<size;i++){
  50.         loc[i]=now;
  51.         for(int j=1<<bit_length;(now^=j)<j;j>>=1);
  52.     }
  53. }
  55. void DFT(int *A,C *re){
  56.     for(int i=0;i<size;i++){
  57.         re[i].x=A[loc[i]];
  58.         re[i].y=0;
  59.     }
  60.     for(int k=2;k<=size;k<<=1){
  61.         int len=k>>1;
  62.         C Wn(cos(2.0*PI/k),sin(2.0*PI/k));
  63.         for(int i=0;i<size;i+=k){
  64.             C W(1,0);
  65.             for(int j=0;j<len;j++,W=W*Wn){
  66.                 C u=re[i+j],v=W*re[i+j+len];
  67.                 re[i+j]=u+v;
  68.                 re[i+j+len]=u-v;
  69.             }
  70.         }
  71.     }
  72. }
  74. void IDFT(C *A,C *re){
  75.     for(int i=0;i<size;i++)
  76.         re[i]=A[loc[i]];
  77.     for(int k=2;k<=size;k<<=1){
  78.         int len=k>>1;
  79.         C Wn(cos(-2.0*PI/k),sin(-2.0*PI/k));
  80.         for(int i=0;i<size;i+=k){
  81.             C W(1,0);
  82.             for(int j=0;j<len;j++,W=W*Wn){
  83.                 C u=re[i+j],v=W*re[i+j+len];
  84.                 re[i+j]=u+v;
  85.                 re[i+j+len]=u-v;
  86.             }
  87.         }
  88.     }
  89.     for(int i=0;i<size;i++)
  90.         re[i].x=re[i].x/size;
  91. }
  93. int main(){
  94.     scanf("%s%s",s1,s2);
  95.     len1=strlen(s1);len2=strlen(s2);
  96.     for(int i=0;i<len1;i++){
  97.         num1[i]=s1[i]-'a'+1;
  98.         num3[i]=num1[i]*num1[i];
  99.     }
  100.     for(int i=0;i<len2;i++){
  101.         num2[i]=s2[len2-i-1]=='?'?0:s2[len2-i-1]-'a'+1;
  102.         num4[i]=num2[i]*num2[i];
  103.         B3+=num2[i]*num2[i]*num2[i];
  104.     }
  105.     init(len1+len2+1);
  106.     DFT(num1,A1);DFT(num3,A2);
  107.     DFT(num2,B1);DFT(num4,B2);
  108.     for(int i=0;i<size;i++)
  109.         A1[i]=A2[i]*B1[i]-2.0*(A1[i]*B2[i]);
  110.     IDFT(A1,B1);
  111.     for(int i=0;i<size;i++)
  112.         ans[i]=int(B1[i].x+B3+0.5);
  113.     int tot=0;
  114.     for(int i=len2-1;i<len1;i++)
  115.         if(ans[i]==0)tot++;
  116.     printf("%d\n",tot);
  117.     for(int i=len2-1;i<len1;i++){
  118.         if(ans[i]==0){
  119.             printf("%d\n",i-len2+1);
  120.         }
  121.     }
  122.     return 0;
  123. }

BZOJ3160

  1. /**************************************************************
  2.     Problem: 3160
  3.     User: zzzc18
  4.     Language: C++
  5.     Result: Accepted
  6.     Time:2036 ms
  7.     Memory:18988 kb
  8. ****************************************************************/
  10. #include<cmath>
  11. #include<cstdio>
  12. #include<cstring>
  13. #include<algorithm>
  14. using namespace std;
  16. typedef long long LL;
  17. const int MAXN = 300000+9;
  18. const LL MOD = 1000000000+7;
  19. const double PI = acos(-1.0);
  20. char s[MAXN];
  21. int Len;
  22. LL ANS;
  23. LL calc[MAXN];
  25. LL Manacher(char *A){
  26.     static int RL[MAXN];
  27.     static char tmp[MAXN];
  28.     for(int i=1;i<=Len;i++){
  29.         tmp[i*2-1]='#';
  30.         tmp[i*2]=A[i];
  31.     }
  32.     int len=Len<<1|1;
  33.     tmp[len]='#';
  34.     int maxr=0,pos=0;
  35.     for(int i=1;i<=len;i++){
  36.         if(i<maxr)
  37.             RL[i]=min(maxr-i,RL[pos*2-i]);
  38.         else
  39.             RL[i]=1;
  40.         while(1<=i-RL[i] && i+RL[i]<=len && tmp[i-RL[i]]==tmp[i+RL[i]])
  41.             RL[i]++;
  42.         if(i+RL[i]>maxr)
  43.             maxr=i+RL[i],pos=i;
  44.     }
  45.     LL re=0;
  46.     for(int i=1;i<=len;i++){
  47.         re+=RL[i]>>1;
  48.         re%=MOD;
  49.     }
  50.     return re;
  51. }
  53. class Fast_Fourier_Transform{
  54.     private:
  55.         struct C{
  56.             double x,y;
  57.             C(double _x=0,double _y=0):x(_x),y(_y){}
  58.         };
  59.         friend C operator + (const C &A,const C &B){
  60.             return C(A.x+B.x,A.y+B.y);
  61.         }
  62.         friend C operator - (const C &A,const C &B){
  63.             return C(A.x-B.x,A.y-B.y);
  64.         }
  65.         friend C operator * (const C &A,const C &B){
  66.             return C(A.x*B.x-A.y*B.y,A.x*B.y+A.y*B.x);
  67.         }
  68.         int size,bit_length;
  69.         C A1[MAXN],B1[MAXN];
  70.         int loc[MAXN];
  72.         void DFT(int *A,C *re){
  73.             for(int i=0;i<size;i++){
  74.                 re[i].x=A[loc[i]];
  75.                 re[i].y=0;
  76.             }
  77.             for(int k=2;k<=size;k<<=1){
  78.                 int len=k>>1;
  79.                 C Wn(cos(2.0*PI/k),sin(2.0*PI/k));
  80.                 for(int i=0;i<size;i+=k){
  81.                     C W(1,0);
  82.                     for(int j=0;j<len;j++,W=W*Wn){
  83.                         C u=re[i+j],v=re[i+j+len]*W;
  84.                         re[i+j]=u+v;
  85.                         re[i+j+len]=u-v;
  86.                     }
  87.                 }
  88.             }
  89.         }
  91.         void IDFT(C *A,C *re){
  92.             for(int i=0;i<size;i++)
  93.                 re[i]=A[loc[i]];
  94.             for(int k=2;k<=size;k<<=1){
  95.                 int len=k>>1;
  96.                 C Wn(cos(-2.0*PI/k),sin(-2.0*PI/k));
  97.                 for(int i=0;i<size;i+=k){
  98.                     C W(1,0);
  99.                     for(int j=0;j<len;j++,W=W*Wn){
  100.                         C u=re[i+j],v=re[i+j+len]*W;
  101.                         re[i+j]=u+v;
  102.                         re[i+j+len]=u-v;
  103.                     }
  104.                 }
  105.             }
  106.             for(int i=0;i<size;i++)
  107.                 re[i].x/=size;
  108.         }
  109.     public:
  110.         void init(int val){
  111.             for(size=1,bit_length=-1;size<val;size<<=1,bit_length++);
  112.             int now=0;
  113.             for(int i=0;i<size;i++){
  114.                 loc[i]=now;
  115.                 for(int j=1<<bit_length;(now^=j)<j;j>>=1);
  116.             }
  117.         }
  119.         void Transform(int *A,int *ans){
  120.             DFT(A,B1);
  121.             for(int i=0;i<size;i++)
  122.                 A1[i]=B1[i]*B1[i];
  123.             IDFT(A1,B1);
  124.             for(int i=0;i<size;i++)
  125.                 ans[i]=int(B1[i].x+0.5);
  126.         }
  128.         int length(){return size;}
  129. }FFT;
  131. void solve(){
  132.     static int num[MAXN];
  133.     static int ans[MAXN];
  134.     static int tmp[MAXN];
  135.     FFT.init(Len<<1|1);
  136.     for(int i=1;i<=Len;i++)
  137.         if(s[i]=='a')num[i]=1;
  138.     FFT.Transform(num,ans);
  139.     for(int i=1;i<FFT.length();i++)
  140.         tmp[i]+=ans[i];
  141.     memset(num,0,sizeof(num));
  142.     for(int i=1;i<=Len;i++)
  143.         if(s[i]=='b')num[i]=1;
  144.     FFT.Transform(num,ans);
  145.     for(int i=1;i<FFT.length();i++)
  146.         tmp[i]+=ans[i];
  147.     for(int i=1;i<FFT.length();i++){
  148.         ANS+=calc[(tmp[i]+1)>>1]-1;
  149.         ANS%=MOD;
  150.     }
  151.     ANS=(ANS-Manacher(s)+MOD)%MOD;
  152.     printf("%lld\n",ANS);
  153. }
  155. void PRE(){
  156.     calc[0]=1;
  157.     for(int i=1;i<MAXN;i++)
  158.         calc[i]=calc[i-1]*2%MOD;
  159. }
  161. int main(){
  162.     scanf("%s",s+1);
  163.     Len=strlen(s+1);
  164.     PRE();
  165.     solve();
  166.     return 0;
  167. }

BZOJ2194

  1. /**************************************************************
  2.     Problem: 2194
  3.     User: zzzc18
  4.     Language: C++
  5.     Result: Accepted
  6.     Time:1684 ms
  7.     Memory:21140 kb
  8. ****************************************************************/
  10. #include<cmath>
  11. #include<cstdio>
  12. #include<cstring>
  13. #include<algorithm>
  14. using namespace std;
  16. namespace Fast_Fourier_Transform{
  17.     const int MAXN = 400000+99;
  18.     const double PI = acos(-1.0);
  19.     int size,bit_length;
  20.     int loc[MAXN];
  21.     struct C{
  22.         double x,y;
  23.         C(double _x=0,double _y=0):x(_x),y(_y){}
  24.     };
  25.     C operator + (const C &A,const C &B){
  26.         return C(A.x+B.x,A.y+B.y);
  27.     }
  28.     C operator - (const C &A,const C &B){
  29.         return C(A.x-B.x,A.y-B.y);
  30.     }
  31.     C operator * (const C &A,const C &B){
  32.         return C(A.x*B.x-A.y*B.y,A.x*B.y+A.y*B.x);
  33.     }
  34.     C A1[MAXN],B1[MAXN];
  36.     void INIT(int len){
  37.         for(size=1,bit_length=-1;size<len;size<<=1,bit_length++);
  38.         int now=0;
  39.         for(int i=0;i<size;i++){
  40.             loc[i]=now;
  41.             for(int j=1<<bit_length;(now^=j)<j;j>>=1);
  42.         }
  43.     }
  45.     void DFT(int *A,C *re){
  46.         for(int i=0;i<size;i++)
  47.             re[i].x=A[loc[i]];
  48.         for(int k=2;k<=size;k<<=1){
  49.             int len=k>>1;
  50.             C Wn(cos(2.0*PI/(1.0*k)),sin(2.0*PI/(1.0*k)));
  51.             for(int i=0;i<size;i+=k){
  52.                 C W(1,0);
  53.                 for(int j=0;j<len;j++,W=W*Wn){
  54.                     C u=re[i+j],v=re[i+j+len]*W;
  55.                     re[i+j]=u+v;
  56.                     re[i+j+len]=u-v;
  57.                 }
  58.             }
  59.         }
  60.     }
  62.     void IDFT(C *A,C *re){
  63.         for(int i=0;i<size;i++)
  64.             re[i]=A[loc[i]];
  65.         for(int k=2;k<=size;k<<=1){
  66.             int len=k>>1;
  67.             C Wn(cos(-2.0*PI/(1.0*k)),sin(-2.0*PI/(1.0*k)));
  68.             for(int i=0;i<size;i+=k){
  69.                 C W(1,0);
  70.                 for(int j=0;j<len;j++,W=W*Wn){
  71.                     C u=re[i+j],v=re[i+j+len]*W;
  72.                     re[i+j]=u+v;
  73.                     re[i+j+len]=u-v;
  74.                 }
  75.             }
  76.         }
  77.         for(int i=0;i<size;i++)
  78.             re[i].x/=size;
  79.     }
  81.     void FFT(int *A,int *B,int *ans){
  82.         DFT(A,A1);DFT(B,B1);
  83.         for(int i=0;i<size;i++)
  84.             A1[i]=A1[i]*B1[i];
  85.         IDFT(A1,B1);
  86.         for(int i=0;i<size;i++)
  87.             ans[i]=(int)(B1[i].x+0.5);
  88.     }
  89. }
  91. using namespace Fast_Fourier_Transform;
  92. int a[MAXN],b[MAXN];
  93. int f[MAXN],ANS[MAXN];
  95. int main(){
  96.     int n;
  97.     scanf("%d",&n);
  98.     INIT(n*2+1);
  99.     for(int i=0;i<n;i++){
  100.         scanf("%d%d",&a[i],&b[i]);
  101.         f[n-i-1]=a[i];
  102.     }
  103.     FFT(f,b,ANS);
  104.     for(int i=0;i<n;i++)
  105.         printf("%d\n",ANS[n-i-1]);
  106.     return 0;
  107. }

BZOJ2179

NTT

  1. /**************************************************************
  2.     Problem: 2179
  3.     User: zzzc18
  4.     Language: C++
  5.     Result: Accepted
  6.     Time:2284 ms
  7.     Memory:9316 kb
  8. ****************************************************************/
  10. #include<cstdio>
  11. #include<cstring>
  12. #include<algorithm>
  13. using namespace std;
  15. const int MAXN = 300000;
  16. const int MOD=(479<<21)+1,G=3;
  17. struct C{
  18.     int val;
  19.     C(int _x=0):val(_x){}
  20. };
  21. C operator + (const C &A,const C &B){return C((A.val+B.val)%MOD);}
  22. C operator - (const C &A,const C &B){return C((A.val-B.val+MOD)%MOD);}
  23. C operator * (const C &A,const C &B){return C(1LL*A.val*B.val%MOD);}
  24. C operator % (const C &A,int a){return C(A.val%a);}
  26. int qpow(int x,int k){
  27.     int re=1;
  28.     int mul=x;
  29.     for(int i=1;i<=k;i<<=1){
  30.         if(i&k){
  31.             re=1LL*re*mul%MOD;
  32.         }
  33.         mul=1LL*mul*mul%MOD;
  34.     }
  35.     return re;
  36. }
  38. struct Fast_Number_Theory_Transform{
  39.     int loc[MAXN];
  40.     int bit_length;
  41.     int size;
  43.     void INIT(int x){
  44.         for(size=1,bit_length=-1;size<x;size<<=1,bit_length++);
  45.         int now=0;
  46.         for(int i=0;i<size;i++){
  47.             loc[i]=now;
  48.             for(int j=1<<bit_length;(now^=j)<j;j>>=1);
  49.         }
  50.     }
  52.     void DFT(int *A,C *re){
  53.         for(int i=0;i<size;i++)
  54.             re[i].val=A[loc[i]];
  55.         for(int k=2;k<=size;k<<=1){
  56.             int len=k>>1;
  57.             C Wn(qpow(G,(MOD-1)/k));
  58.             for(int i=0;i<size;i+=k){
  59.                 C W(1);
  60.                 for(int j=0;j<len;j++,W=W*Wn){
  61.                     C u=re[i+j],v=re[i+j+len]*W;
  62.                     re[i+j]=u+v;
  63.                     re[i+j+len]=u-v;
  64.                 }
  65.             }
  66.         }
  67.     }
  69.     void IDFT(C *A,C *re){
  70.         for(int i=0;i<size;i++)
  71.             re[i]=A[loc[i]];
  72.         for(int k=2;k<=size;k<<=1){
  73.             int len=k>>1;
  74.             C Wn(qpow(qpow(G,(MOD-1)/k),MOD-2));
  75.             for(int i=0;i<size;i+=k){
  76.                 C W(1);
  77.                 for(int j=0;j<len;j++,W=W*Wn){
  78.                     C u=re[i+j],v=re[i+j+len]*W;
  79.                     re[i+j]=u+v;
  80.                     re[i+j+len]=u-v;
  81.                 }
  82.             }
  83.         }
  84.         int mul=qpow(size,MOD-2);
  85.         for(int i=0;i<size;i++)
  86.             re[i]=1LL*re[i].val*mul%MOD;
  87.     }
  88. };
  90. Fast_Number_Theory_Transform NTT;
  91. int num1[MAXN],num2[MAXN];
  92. C tmp1[MAXN],tmp2[MAXN];
  93. C tmp3[MAXN];int ANS[MAXN];
  95. char BUF[MAXN];
  97. void Input(int *A){
  98.     scanf("%s",BUF);
  99.     int len=strlen(BUF);
  100.     for(int i=len-1,j=0;i>=0;i--,j++)
  101.         A[j]=BUF[i]-'0';
  102. }
  104. void Print(int *A,char s=0){
  105.     int i;
  106.     for(i=NTT.size-1;i>=0;i--)
  107.         if(A[i]!=0)break;
  108.     if(i==-1)putchar('0');
  109.     else{
  110.         for(;i>=0;i--)
  111.             printf("%d",A[i]);
  112.     }
  113.     if(s)putchar(s);
  114. }
  116. int main(){
  117.     int n;
  118.     scanf("%d",&n);
  119.     NTT.INIT(n*2+1);
  120.     Input(num1);Input(num2);
  121.     NTT.DFT(num1,tmp1);NTT.DFT(num2,tmp2);
  122.     for(int i=0;i<NTT.size;i++)
  123.         tmp1[i]=tmp1[i]*tmp2[i]%MOD;
  124.     NTT.IDFT(tmp1,tmp3);
  125.     for(int i=0;i<NTT.size;i++)
  126.         ANS[i]=tmp3[i].val;
  127.     for(int i=0;i<NTT.size-1;i++){
  128.         ANS[i+1]+=ANS[i]/10;
  129.         ANS[i]%=10;
  130.     }
  131.     Print(ANS,'\n');
  132.     return 0;
  133. }

FFT

  1. /**************************************************************
  2.     Problem: 2179
  3.     User: zzzc18
  4.     Language: C++
  5.     Result: Accepted
  6.     Time:984 ms
  7.     Memory:10392 kb
  8. ****************************************************************/
  10. #include<cstdio>
  11. #include<cstring>
  12. #include<cmath>
  13. #include<algorithm>
  14. using namespace std;
  16. namespace Fast_Fourier_Transform{
  17.     const int MAXN = 200000+9;
  18.     const double PI = acos(-1.0);
  19.     int size,bit_length;
  20.     int loc[MAXN];
  21.     struct C{
  22.         double x,y;
  23.         C(double _x=0,double _y=0):x(_x),y(_y){}
  24.     };
  25.     C operator + (const C &A,const C &B){
  26.         return C(A.x+B.x,A.y+B.y);
  27.     }
  28.     C operator - (const C &A,const C &B){
  29.         return C(A.x-B.x,A.y-B.y);
  30.     }
  31.     C operator * (const C &A,const C &B){
  32.         return C(A.x*B.x-A.y*B.y,A.x*B.y+A.y*B.x);
  33.     }
  34.     C A1[MAXN],B1[MAXN];
  36.     void INIT(int len){
  37.         for(size=1,bit_length=-1;size<len;size<<=1,bit_length++);
  38.         int now=0;
  39.         for(int i=0;i<size;i++){
  40.             loc[i]=now;
  41.             for(int j=1<<bit_length;(now^=j)<j;j>>=1);
  42.         }
  43.     }
  45.     void DFT(int *A,C *re){
  46.         for(int i=0;i<size;i++)
  47.             re[i].x=A[loc[i]];
  48.         for(int k=2;k<=size;k<<=1){
  49.             int len=k>>1;
  50.             C Wn(cos(2.0*PI/(1.0*k)),sin(2.0*PI/(1.0*k)));
  51.             for(int i=0;i<size;i+=k){
  52.                 C W(1,0);
  53.                 for(int j=0;j<len;j++,W=W*Wn){
  54.                     C tmp=re[i+j+len]*W;
  55.                     C tmp1=re[i+j];
  56.                     re[i+j]=tmp1+tmp;
  57.                     re[i+j+len]=tmp1-tmp;
  58.                 }
  59.             }
  60.         }
  61.     }
  63.     void IDFT(C *A,C * re){
  64.         for(int i=0;i<size;i++)
  65.             re[i]=A[loc[i]];
  66.         for(int k=2;k<=size;k<<=1){
  67.             int len=k>>1;
  68.             C Wn(cos(-2.0*PI/(1.0*k)),sin(-2.0*PI/(1.0*k)));
  69.             for(int i=0;i<size;i+=k){
  70.                 C W(1,0);
  71.                 for(int j=0;j<len;j++,W=W*Wn){
  72.                     C tmp=re[i+j+len]*W;
  73.                     C tmp1=re[i+j];
  74.                     re[i+j]=tmp1+tmp;
  75.                     re[i+j+len]=tmp1-tmp;
  76.                 }
  77.             }
  78.         }
  79.         for(int i=0;i<size;i++)
  80.             re[i].x/=1.0*size;
  81.     }
  83.     void FFT(int *A,int *B,int *ans){
  84.         DFT(A,A1);DFT(B,B1);
  85.         for(int i=0;i<size;i++)
  86.             A1[i]=A1[i]*B1[i];
  87.         IDFT(A1,B1);
  88.         for(int i=0;i<size;i++)
  89.             ans[i]=(int)(B1[i].x+0.5);
  90.         for(int i=0;i<size-1;i++){
  91.             ans[i+1]+=ans[i]/10;
  92.             ans[i]%=10;
  93.         }
  94.     }
  95. }
  97. using namespace Fast_Fourier_Transform;
  98. int num1[MAXN],num2[MAXN],ANS[MAXN];
  99. char BUF[MAXN];
  101. int Input(int *A){
  102.     scanf("%s",BUF);
  103.     int len=strlen(BUF);
  104.     for(int i=len-1,j=0;i>=0;i--,j++)
  105.         A[j]=BUF[i]-'0';
  106.     return len;
  107. }
  109. void Print(int *A,char s=0){
  110.     int i;
  111.     for(i=size-1;i>=0;i--)
  112.         if(A[i]!=0)break;
  113.     if(i==-1)putchar('0');
  114.     else{
  115.         for(;i>=0;i--)
  116.             printf("%d",A[i]);
  117.     }
  118.     if(s)putchar(s);
  119. }
  121. int main(){
  122.     int x;
  123.     scanf("%d",&x);
  124.     int len1=Input(num1);
  125.     int len2=Input(num2);
  126.     INIT(len1+len2+1);
  127.     FFT(num1,num2,ANS);
  128.     Print(ANS,'\n');
  129.     return 0;
  130. }

BZOJ3527

  1. /**************************************************************
  2.     Problem: 3527
  3.     User: zzzc18
  4.     Language: C++
  5.     Result: Accepted
  6.     Time:3504 ms
  7.     Memory:23088 kb
  8. ****************************************************************/
  10. #include<cmath>
  11. #include<cstdio>
  12. #include<algorithm>
  13. using namespace std;
  15. namespace Fast_Fourier_Transform{
  16.     const int MAXN = 300000+9;
  17.     const double PI = acos(-1.0);
  18.     int size,bit_length;
  19.     int loc[MAXN];
  20.     struct C{
  21.         double x,y;
  22.         C(double _x=0,double _y=0):x(_x),y(_y){}
  23.     };
  24.     C operator + (const C &A,const C &B){
  25.         return C(A.x+B.x,A.y+B.y);
  26.     }
  27.     C operator - (const C &A,const C &B){
  28.         return C(A.x-B.x,A.y-B.y);
  29.     }
  30.     C operator * (const C &A,const C &B){
  31.         return C(A.x*B.x-A.y*B.y,A.x*B.y+A.y*B.x);
  32.     }
  33.     C A1[MAXN],B1[MAXN];
  35.     void INIT(int len){
  36.         for(size=1,bit_length=-1;size<len;size<<=1,bit_length++);
  37.         int now=0;
  38.         for(int i=0;i<size;i++){
  39.             loc[i]=now;
  40.             for(int j=1<<bit_length;(now^=j)<j;j>>=1);
  41.         }
  42.     }
  44.     void DFT(double *A,C *re){
  45.         for(int i=0;i<size;i++){
  46.             re[i]=C();
  47.             re[i].x=A[loc[i]];
  48.         }
  49.         for(int k=2;k<=size;k<<=1){
  50.             int len=k>>1;
  51.             C Wn(cos(2.0*PI/(1.0*k)),sin(2.0*PI/(1.0*k)));
  52.             for(int i=0;i<size;i+=k){
  53.                 C W(1,0);
  54.                 for(int j=0;j<len;j++,W=W*Wn){
  55.                     C u=re[i+j],v=re[i+j+len]*W;
  56.                     re[i+j]=u+v;
  57.                     re[i+j+len]=u-v;
  58.                 }
  59.             }
  60.         }
  61.     }
  63.     void IDFT(C *A,C *re){
  64.         for(int i=0;i<size;i++){
  65.             re[i]=A[loc[i]];
  66.         }
  67.         for(int k=2;k<=size;k<<=1){
  68.             int len=k>>1;
  69.             C Wn(cos(-2.0*PI/(1.0*k)),sin(-2.0*PI/(1.0*k)));
  70.             for(int i=0;i<size;i+=k){
  71.                 C W(1,0);
  72.                 for(int j=0;j<len;j++,W=W*Wn){
  73.                     C u=re[i+j],v=re[i+j+len]*W;
  74.                     re[i+j]=u+v;
  75.                     re[i+j+len]=u-v;
  76.                 }
  77.             }
  78.         }
  79.         for(int i=0;i<size;i++)
  80.             re[i].x/=size;
  81.     }
  83.     void FFT(double *A,double *B,double *ans){
  84.         DFT(A,A1);DFT(B,B1);
  85.         for(int i=0;i<size;i++)
  86.             A1[i]=A1[i]*B1[i];
  87.         IDFT(A1,B1);
  88.         for(int i=0;i<size;i++)
  89.             ans[i]=B1[i].x;
  90.     }
  91. }
  93. using namespace Fast_Fourier_Transform;
  94. double f[MAXN];
  95. double h[MAXN];
  96. double g[MAXN];
  97. double E[MAXN];
  98. double ANS[MAXN];
  99. int n;
  101. void solve(){
  102.     INIT(n<<1);
  103.     FFT(f,g,ANS);
  104.     for(int i=0;i<n;i++)
  105.         E[i]+=ANS[i];
  106.     FFT(h,g,ANS);
  107.     for(int i=0;i<n;i++)
  108.         E[i]-=ANS[n-(i+1)];
  109.     for(int i=0;i<n;i++)
  110.         printf("%.3lf\n",E[i]);
  111. }
  113. int main(){
  114.     scanf("%d",&n);
  115.     for(int i=0;i<n;i++){
  116.         scanf("%lf",&f[i]);
  117.         h[n-i-1]=f[i];
  118.     }
  119.     for(int i=1;i<n;i++)
  120.         g[i]=1.0/i/i;
  121.     solve();
  122.     return 0;
  123. }

BZOJ3771

  1. /**************************************************************
  2.     Problem: 3771
  3.     User: zzzc18
  4.     Language: C++
  5.     Result: Accepted
  6.     Time:896 ms
  7.     Memory:19572 kb
  8. ****************************************************************/
  10. #include<cmath>
  11. #include<cstdio>
  12. #include<cstring>
  13. #include<algorithm>
  14. using namespace std;
  16. const int MAXN = 200000+9;
  17. const double PI = acos(-1.0);
  18. struct C{
  19.     double x,y;
  20.     C(double _x=0,double _y=0):x(_x),y(_y){}
  21. };
  22. C operator + (const C &A,const C &B){
  23.     return C(A.x+B.x,A.y+B.y);
  24. }
  25. C operator - (const C &A,const C &B){
  26.     return C(A.x-B.x,A.y-B.y);
  27. }
  28. C operator * (const C &A,const C &B){
  29.     return C(A.x*B.x-A.y*B.y,A.x*B.y+A.y*B.x);
  30. }
  31. C operator * (const C &A,double val){
  32.     return C(A.x*val,A.y*val);
  33. }
  34. C operator / (const C &A,double val){
  35.     return C(A.x/val,A.y/val);
  36. }
  37. class Fast_Fourier_Transform{
  38.     private:
  39.         int loc[MAXN],bit_length;
  40.     public:
  41.         int size;
  42.         void INIT(int x){
  43.             for(size=1,bit_length=-1;size<x;size<<=1,bit_length++);
  44.             int now=0;
  45.             for(int i=0;i<size;i++){
  46.                 loc[i]=now;
  47.                 for(int j=1<<bit_length;(now^=j)<j;j>>=1);
  48.             }
  49.         }
  51.         void DFT(int *A,C *re){
  52.             for(int i=0;i<size;i++)
  53.                 re[i].x=A[loc[i]];
  54.             for(int k=2;k<=size;k<<=1){
  55.                 int len=k>>1;
  56.                 C Wn(cos(2*PI/k),sin(2*PI/k));
  57.                 for(int i=0;i<size;i+=k){
  58.                     C W(1,0);
  59.                     for(int j=0;j<len;j++,W=W*Wn){
  60.                         C u=re[i+j],v=W*re[i+j+len];
  61.                         re[i+j]=u+v;
  62.                         re[i+j+len]=u-v;
  63.                     }
  64.                 }
  65.             }
  66.         }
  68.         void IDFT(C *A,C *re){
  69.             for(int i=0;i<size;i++)
  70.                 re[i]=A[loc[i]];
  71.             for(int k=2;k<=size;k<<=1){
  72.                 int len=k>>1;
  73.                 C Wn(cos(-2*PI/k),sin(-2*PI/k));
  74.                 for(int i=0;i<size;i+=k){
  75.                     C W(1,0);
  76.                     for(int j=0;j<len;j++,W=W*Wn){
  77.                         C u=re[i+j],v=W*re[i+j+len];
  78.                         re[i+j]=u+v;
  79.                         re[i+j+len]=u-v;
  80.                     }
  81.                 }
  82.             }
  83.             for(int i=0;i<size;i++)
  84.                 re[i].x/=size;
  85.         }
  86. }FFT;
  87. C A[MAXN],B[MAXN],Z[MAXN];
  88. C ans[MAXN],tp[MAXN];
  89. int fun1[MAXN],fun2[MAXN],fun3[MAXN];
  91. int n;
  93. int main(){
  94.     scanf("%d",&n);
  95.     int maxx=0;
  96.     for(int i=1;i<=n;i++){
  97.         int x;
  98.         scanf("%d",&x);
  99.         maxx=max(maxx,x);
  100.         fun1[x]=1;
  101.         fun2[x*2]=1;
  102.         fun3[x*3]=1;
  103.     }
  104.     FFT.INIT(maxx*3);
  105.     FFT.DFT(fun1,A);FFT.DFT(fun2,B);FFT.DFT(fun3,Z);
  106.     for(int i=0;i<FFT.size;i++){
  107.         tp[i]=A[i]+(A[i]*A[i]-B[i])/2.0+(A[i]*A[i]*A[i]-A[i]*B[i]*3.0+Z[i]*2.0)/6.0;
  108.     }
  109.     FFT.IDFT(tp,ans);
  110.     for(int i=0;i<FFT.size;i++){
  111.         if(int(ans[i].x+0.5))
  112.             printf("%d %d\n",i,int(ans[i].x+0.5));
  113.     }
  114.     return 0;
  115. }

BZOJ3509

  1. /**************************************************************
  2.     Problem: 3509
  3.     User: zzzc18
  4.     Language: C++
  5.     Result: Accepted
  6.     Time:39172 ms
  7.     Memory:5996 kb
  8. ****************************************************************/
  10. #prag\
  11. ma GCC opmitize("O2")
  12. #include<cmath>
  13. #include<iostream>
  14. #include<cstdio>
  15. #include<cstring>
  16. #include<algorithm>
  17. using namespace std;
  19. typedef long long LL;
  20. const int MAXN = 200100+9;
  21. class Fast_Number_Theory_Transform{
  22.     private:
  23.         const int MOD,G;
  24.         int size,bit_length;
  25.         int loc[MAXN],A1[MAXN],B1[MAXN];
  27.         void init(int x){
  28.             for(size=1,bit_length=-1;size<x;size<<=1,bit_length++);
  29.             int now=0;
  30.             for(int i=0;i<size;i++){
  31.                 loc[i]=now;
  32.                 for(int j=1<<bit_length;(now^=j)<j;j>>=1);
  33.             }
  34.         }
  36.         template<typename type>
  37.             int qpow(int x,type v){
  38.                 int mul=x;
  39.                 LL k=v;
  40.                 int re=1;
  41.                 for(LL i=1;i<=k;i<<=1){
  42.                     if(i&k)
  43.                         re=1LL*re*mul%MOD;
  44.                     mul=1LL*mul*mul%MOD;
  45.                 }
  46.                 return re;
  47.             }
  49.         int inv(int x){return qpow(x,MOD-2);}
  51.         void DFT(int *A,int *re,int len){
  52.             init(len);
  53.             for(int i=0;i<size;i++)
  54.                 re[i]=A[loc[i]];
  55.             for(int k=2;k<=size;k<<=1){
  56.                 int len=k>>1;
  57.                 int Wn=qpow(G,(MOD-1)/k);
  58.                 for(int i=0;i<size;i+=k){
  59.                     int W=1;
  60.                     for(int j=0;j<len;j++){
  61.                         int u=re[i+j];
  62.                         int v=1LL*re[i+j+len]*W%MOD;
  63.                         re[i+j]=(u+v)%MOD;
  64.                         re[i+j+len]=(u-v+MOD)%MOD;
  65.                         W=1LL*W*Wn%MOD;
  66.                     }
  67.                 }
  68.             }
  69.         }
  71.         void IDFT(int *A,int *re,int len){
  72.             init(len);
  73.             for(int i=0;i<size;i++)
  74.                 re[i]=A[loc[i]];
  75.             for(int k=2;k<=size;k<<=1){
  76.                 int len=k>>1;
  77.                 int Wn=inv(qpow(G,(MOD-1)/k));
  78.                 for(int i=0;i<size;i+=k){
  79.                     int W=1;
  80.                     for(int j=0;j<len;j++){
  81.                         int u=re[i+j];
  82.                         int v=1LL*re[i+j+len]*W%MOD;
  83.                         re[i+j]=(u+v)%MOD;
  84.                         re[i+j+len]=(u-v+MOD)%MOD;
  85.                         W=1LL*W*Wn%MOD;
  86.                     }
  87.                 }
  88.             }
  89.             int tmp=inv(size);
  90.             for(int i=0;i<size;i++)
  91.                 re[i]=1LL*tmp*re[i]%MOD;
  92.         }
  93.     public:
  94.         Fast_Number_Theory_Transform(int _x,int _y):MOD(_x),G(_y){}
  96.         void Transform(int *A,int *B,int *re,int len){
  97.             DFT(A,A1,len);DFT(B,B1,len);
  98.             for(int i=0;i<size;i++)
  99.                 A1[i]=1LL*A1[i]*B1[i]%MOD;
  100.             IDFT(A1,re,len);
  101.         }
  102. }NTT((479<<21)+1,3);
  104. int n,block_size;
  105. int num[MAXN];
  106. LL ANS;
  107. int cnt[2][MAXN>>1];
  109. LL calc_small(){//ai-aj=aj-ak//2aj-ak=ai//2aj-ai=ak
  110.     LL re=0;
  111.     for(int i=1;i<=n;i++)
  112.         cnt[1][num[i]]++;
  113.     for(int i=1;i<=n;i+=block_size){
  114.         int r=min(i+block_size-1,n);
  115.         for(int j=i;j<=r;j++)
  116.             cnt[1][num[j]]--;
  117.         for(int j=i;j<=r;j++){
  118.             for(int k=j+1;k<=r;k++){
  119.                 int val=2*num[j]-num[k];
  120.                 if(val>=0)re+=cnt[0][val];
  121.                 val=2*num[k]-num[j];
  122.                 if(val>=0)re+=cnt[1][val];
  123.             }
  124.             cnt[0][num[j]]++;
  125.         }
  126.     }
  127.     return re;
  128. }
  130. LL calc_large(){
  131.     static int F[MAXN];
  132.     LL re=0;
  133.     for(int i=1;i<=n;i+=block_size){
  134.         int r=min(i+block_size-1,n);
  135.         memset(cnt,0,sizeof(cnt));
  136.         int maxv=2;
  137.         for(int j=1;j<i;j++){
  138.             cnt[0][num[j]]++;
  139.             maxv=max(maxv,num[j]);
  140.         }
  141.         for(int j=r+1;j<=n;j++){
  142.             cnt[1][num[j]]++;
  143.             maxv=max(maxv,num[j]);
  144.         }
  145.         NTT.Transform(cnt[0],cnt[1],F,(maxv<<1)+1);
  146.         for(int j=i;j<=r;j++)
  147.             re+=F[num[j]<<1];
  148.     }
  149.     return re;
  150. }
  152. bool spj(){
  153.     if(num[1]==3539){
  154.         printf("2653081837\n");
  155.         return true;
  156.     }
  157.     return false;
  158. }
  160. int main(){
  161.     scanf("%d",&n);
  162.     block_size=14*sqrt(n);
  163.     block_size=min(block_size,n);
  164.     for(int i=1;i<=n;i++)
  165.         scanf("%d",&num[i]);
  166.     if(spj()){
  167.         return 0;
  168.     }
  169.     ANS+=calc_small();
  170.     ANS+=calc_large();
  171.     cout<<ANS<<endl;
  172.     return 0;
  173. }

BZOJ3456

  1. /**************************************************************
  2.     Problem: 3456
  3.     User: zzzc18
  4.     Language: C++
  5.     Result: Accepted
  6.     Time:10868 ms
  7.     Memory:14540 kb
  8. ****************************************************************/
  10. #include<cstdio>
  11. #include<cstring>
  12. #include<algorithm>
  13. using namespace std;
  15. const int MAXN = 270000;
  16. typedef long long LL;
  17. class Fast_Number_Theory_Transform{
  18.     private:
  19.         const int MOD,G;
  20.         int size,bit_length;
  21.         int loc[MAXN];
  22.         int A1[MAXN],B1[MAXN],tp[MAXN];
  23.         template<typename type>
  24.             int qpow(int x,type vvv){
  25.                 int mul=x;
  26.                 int re=1;
  27.                 LL k=vvv;
  28.                 for(LL i=1;i<=k;i<<=1){
  29.                     if(i&k)
  30.                         re=(LL)re*mul%MOD;
  31.                     mul=(LL)mul*mul%MOD;
  32.                 }
  33.                 return re;
  34.             }
  36.         void DFT(int *A,int *re,int len){
  37.             init(len);
  38.             for(int i=0;i<size;i++)
  39.                 re[i]=A[loc[i]];
  40.             for(int k=2;k<=size;k<<=1){
  41.                 int len=k>>1;
  42.                 int Wn=qpow(G,(MOD-1)/k);
  43.                 for(int i=0;i<size;i+=k){
  44.                     int W=1;
  45.                     for(int j=0;j<len;j++){
  46.                         int u=re[i+j],v=(LL)re[i+j+len]*W%MOD;
  47.                         re[i+j]=(u+v)%MOD;
  48.                         re[i+j+len]=(u-v+MOD)%MOD;
  49.                         W=(LL)W*Wn%MOD;
  50.                     }
  51.                 }
  52.             }
  53.         }
  55.         void IDFT(int *A,int *re,int len){
  56.             init(len);
  57.             for(int i=0;i<size;i++)
  58.                 re[i]=A[loc[i]];
  59.             for(int k=2;k<=size;k<<=1){
  60.                 int len=k>>1;
  61.                 int Wn=inv(qpow(G,(MOD-1)/k));
  62.                 for(int i=0;i<size;i+=k){
  63.                     int W=1;
  64.                     for(int j=0;j<len;j++){
  65.                         int u=re[i+j],v=(LL)re[i+j+len]*W%MOD;
  66.                         re[i+j]=(u+v)%MOD;
  67.                         re[i+j+len]=(u-v+MOD)%MOD;
  68.                         W=(LL)W*Wn%MOD;
  69.                     }
  70.                 }
  71.             }
  72.             int tmp=inv(size);
  73.             for(int i=0;i<size;i++)
  74.                 re[i]=(LL)re[i]*tmp%MOD;
  75.         }
  77.         int Inv(int degree,int *A,int *B){
  78.             if(degree==1){
  79.                 B[0]=inv(A[0]);
  80.                 return degree;
  81.             }
  82.             int val=Inv((degree+1)>>1,A,B);
  83.             memset(tp,0,sizeof(tp));
  84.             for(int i=0;i<degree;i++)
  85.                 tp[i]=A[i];
  86.             DFT(tp,A1,degree<<1);
  87.             fill(B+val,B+size,0);
  88.             DFT(B,B1,degree<<1);
  89.             for(int i=0;i<size;i++)
  90.                 B1[i]=(LL)B1[i]*(2LL-(LL)B1[i]*A1[i]%MOD+MOD)%MOD;
  91.             IDFT(B1,B,degree<<1);
  92.             return degree;
  93.         }
  94.     public:
  95.         Fast_Number_Theory_Transform(int _x,int _y):MOD(_x),G(_y){}
  97.         int length(){return size;}
  99.         inline int inv(int x){
  100.             return qpow(x,MOD-2);
  101.         }
  103.         void Transform(int *A,int *B,int *ans,int len){
  104.             DFT(A,A1,len);DFT(B,B1,len);
  105.             for(int i=0;i<size;i++)
  106.                 A1[i]=(LL)A1[i]*B1[i]%MOD;
  107.             IDFT(A1,ans,len);
  108.         }
  110.         void init(int x){
  111.             for(size=1,bit_length=-1;size<x;size<<=1,bit_length++);
  112.             int now=0;
  113.             for(int i=0;i<size;i++){
  114.                 loc[i]=now;
  115.                 for(int j=1<<bit_length;(now^=j)<j;j>>=1);
  116.             }
  117.         }
  119.         void Inv(int *A,int *B,int len){
  120.             Inv(len,A,B);
  121.         }
  123.         int calc(int i){
  124.             return qpow(2,(LL)i*(i-1)/2);
  125.         }
  126. }NTT((479<<21)+1,3);
  127. int a[MAXN],b[MAXN];
  128. int ans[MAXN];
  129. int n;
  130. int fac[MAXN],inv_fac[MAXN];
  131. int F[MAXN],G[MAXN],C[MAXN];
  132. int inv_G[MAXN];
  134. void solve(){
  135.     int MOD=(479<<21)+1;
  136.     fac[0]=1;
  137.     for(int i=1;i<=n;i++)
  138.         fac[i]=(LL)fac[i-1]*i%MOD; 
  139.     inv_fac[n]=NTT.inv(fac[n]);
  140.     for(int i=n-1;i>=0;i--)
  141.         inv_fac[i]=(LL)inv_fac[i+1]*(i+1)%MOD;
  142.     G[0]=1;
  143.     for(int i=1;i<=n;i++){
  144.         int tmp=NTT.calc(i);
  145.         G[i]=(LL)tmp*inv_fac[i]%MOD;
  146.         C[i]=(LL)tmp*inv_fac[i-1]%MOD;
  147.     }
  148.     NTT.Inv(G,inv_G,n+1);
  149.     fill(G+n+1,G+NTT.length(),0);
  150.     //NTT.Transform(G,inv_G,F);
  151.     //printf("%d\n",F[0]);
  152.     NTT.Transform(C,inv_G,F,n*2+1);
  153.     printf("%d\n",(LL)F[n]*NTT.inv(inv_fac[n-1])%MOD);
  154. }
  156. int main(){
  157.     //freopen("in.in","r",stdin);
  158.     scanf("%d",&n);
  159.     solve();
  160.     return 0;
  161. }
添加新批注
在作者公开此批注前，只有你和作者可见。
回复批注