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

快速傅里叶变换学习笔记

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. }
  13. int FindG(int val){
  14. int re=1;
  15. while(!judge(re)) re++;
  16. return re;
  17. }

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

就变成了

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

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


的第 位就是最终答案

时间复杂度

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

BZOJ3992

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

BZOJ4503

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

BZOJ3160

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

BZOJ2194

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

BZOJ2179

NTT

  1. /**************************************************************
  2. Problem: 2179
  3. User: zzzc18
  4. Language: C++
  5. Result: Accepted
  6. Time:2284 ms
  7. Memory:9316 kb
  8. ****************************************************************/
  9. #include<cstdio>
  10. #include<cstring>
  11. #include<algorithm>
  12. using namespace std;
  13. const int MAXN = 300000;
  14. const int MOD=(479<<21)+1,G=3;
  15. struct C{
  16. int val;
  17. C(int _x=0):val(_x){}
  18. };
  19. C operator + (const C &A,const C &B){return C((A.val+B.val)%MOD);}
  20. C operator - (const C &A,const C &B){return C((A.val-B.val+MOD)%MOD);}
  21. C operator * (const C &A,const C &B){return C(1LL*A.val*B.val%MOD);}
  22. C operator % (const C &A,int a){return C(A.val%a);}
  23. int qpow(int x,int k){
  24. int re=1;
  25. int mul=x;
  26. for(int i=1;i<=k;i<<=1){
  27. if(i&k){
  28. re=1LL*re*mul%MOD;
  29. }
  30. mul=1LL*mul*mul%MOD;
  31. }
  32. return re;
  33. }
  34. struct Fast_Number_Theory_Transform{
  35. int loc[MAXN];
  36. int bit_length;
  37. int size;
  38. void INIT(int x){
  39. for(size=1,bit_length=-1;size<x;size<<=1,bit_length++);
  40. int now=0;
  41. for(int i=0;i<size;i++){
  42. loc[i]=now;
  43. for(int j=1<<bit_length;(now^=j)<j;j>>=1);
  44. }
  45. }
  46. void DFT(int *A,C *re){
  47. for(int i=0;i<size;i++)
  48. re[i].val=A[loc[i]];
  49. for(int k=2;k<=size;k<<=1){
  50. int len=k>>1;
  51. C Wn(qpow(G,(MOD-1)/k));
  52. for(int i=0;i<size;i+=k){
  53. C W(1);
  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. }
  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(qpow(qpow(G,(MOD-1)/k),MOD-2));
  68. for(int i=0;i<size;i+=k){
  69. C W(1);
  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. int mul=qpow(size,MOD-2);
  78. for(int i=0;i<size;i++)
  79. re[i]=1LL*re[i].val*mul%MOD;
  80. }
  81. };
  82. Fast_Number_Theory_Transform NTT;
  83. int num1[MAXN],num2[MAXN];
  84. C tmp1[MAXN],tmp2[MAXN];
  85. C tmp3[MAXN];int ANS[MAXN];
  86. char BUF[MAXN];
  87. void Input(int *A){
  88. scanf("%s",BUF);
  89. int len=strlen(BUF);
  90. for(int i=len-1,j=0;i>=0;i--,j++)
  91. A[j]=BUF[i]-'0';
  92. }
  93. void Print(int *A,char s=0){
  94. int i;
  95. for(i=NTT.size-1;i>=0;i--)
  96. if(A[i]!=0)break;
  97. if(i==-1)putchar('0');
  98. else{
  99. for(;i>=0;i--)
  100. printf("%d",A[i]);
  101. }
  102. if(s)putchar(s);
  103. }
  104. int main(){
  105. int n;
  106. scanf("%d",&n);
  107. NTT.INIT(n*2+1);
  108. Input(num1);Input(num2);
  109. NTT.DFT(num1,tmp1);NTT.DFT(num2,tmp2);
  110. for(int i=0;i<NTT.size;i++)
  111. tmp1[i]=tmp1[i]*tmp2[i]%MOD;
  112. NTT.IDFT(tmp1,tmp3);
  113. for(int i=0;i<NTT.size;i++)
  114. ANS[i]=tmp3[i].val;
  115. for(int i=0;i<NTT.size-1;i++){
  116. ANS[i+1]+=ANS[i]/10;
  117. ANS[i]%=10;
  118. }
  119. Print(ANS,'\n');
  120. return 0;
  121. }

FFT

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

BZOJ3527

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

BZOJ3771

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

BZOJ3509

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

BZOJ3456

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