干货满满的良心博客
- 传送至只有小怪的村庄——请开始你的逆天之路
- A:P1919
- B:P4157
- 刷怪升级——转战玄灵大陆
- C:P6300
- D:P3763
- E:P3321
- F:P5641
- G:P4986
- H:P4721——获得知识药剂一瓶——分治FFTFFTFFT
- I:P3338
- J:P4173
- K:P5748
- L:P3702
- M:P5667
- N:P4245——拾得知识药丸一颗——三模数NTTNTTNTT
- O:P5488
- P:P4199
- Q:P4061
- 缘遇逆天秘籍——金牌打野
- R:P4841
- S:[UVA4671]K-neighbor substrings
- T:SP8372 TSUM - Triple Sums
- U:[CodeChef - COUNTARI]Arithmetic Progressions
- V:CodeForces528D Fuzzy Search
- W:CodeForces 954I Yet Another String Matching Problem
因为对卷积还不是很清楚,所以开始疯狂的刷题之路,会慢慢更新的。
前面没有来源标注的题目题号洛谷上面都有。
FFT,NTTFFT,NTTFFT,NTT是卷积运算中常见的优化→O(nlogn)\rightarrow O(nlogn)→O(nlogn)
因为文字量极大,很有可能会有细节错误,欢迎指出,谢谢o(=·ω·=)m
传送至只有小怪的村庄——请开始你的逆天之路
A:P1919
【模板】A*B Problem升级版(FFT快速傅里叶)
将输入的a,ba,ba,b每一位拆成对应的多项式系数
手玩一下普通乘法的计算法则,发现从左到右对应多项式次数从0,n−10,n-10,n−1依次递增
两数相乘是符合卷积形式,下标是相加的
所以将这个串反转,从左到右从000开始编号
然后套FFTFFTFFT跑出乘积结果,再转回系数表达式
此时某些位置上的值可能很大,这就涉及到进位的问题
从左到右扫一遍即可完成进位
最后倒着输出便是ACACAC
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
#define maxn 4000100
struct complex {double x, i;complex(){}complex( double X, double I ) {x = X, i = I;}
}A[maxn], B[maxn];double pi = acos( -1.0 );complex operator + ( complex a, complex b ) {return complex( a.x + b.x, a.i + b.i );
}complex operator - ( complex a, complex b ) {return complex( a.x - b.x, a.i - b.i );
}complex operator * ( complex a, complex b ) {return complex( a.x * b.x - a.i * b.i, a.x * b.i + a.i * b.x );
}char a[maxn], b[maxn];
int r[maxn], ans[maxn];
int len;void FFT( complex *v, int opt ) {for( int i = 0;i < len;i ++ )if( i < r[i] ) swap( v[i], v[r[i]] );for( int i = 1;i < len;i <<= 1 ) {complex omega( cos( pi / i ), opt * sin( pi / i ) );for( int j = 0;j < len;j += ( i << 1 ) ) {complex w( 1, 0 );for( int k = 0;k < i;k ++, w = w * omega ) {complex x = v[j + k], y = v[j + k + i] * w;v[j + k] = x + y;v[j + k + i] = x - y;}}}
}int main() {scanf( "%s %s", a, b );int n = strlen( a ), m = strlen( b );for( int i = 0;i < n;i ++ ) A[n - i - 1] = complex( a[i] - '0', 0 );for( int i = 0;i < m;i ++ ) B[m - i - 1] = complex( b[i] - '0', 0 );len = 1; int l = 0;while( len < ( n + m ) ) len <<= 1, l ++;for( int i = 0;i < len;i ++ )r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );FFT( A, 1 );FFT( B, 1 );for( int i = 0;i < len;i ++ )A[i] = A[i] * B[i]; FFT( A, -1 );for( int i = 0;i < len;i ++ )ans[i] = ( A[i].x / len + 0.5 );for( int i = 0;i < len;i ++ )if( ans[i] > 9 ) ans[i + 1] += ans[i] / 10, ans[i] %= 10;while( ! ans[len - 1] ) len --;for( int i = len - 1;~ i;i -- )printf( "%d", ans[i] );return 0;
}
B:P4157
[SCOI2006]整数划分
这道题显然大整数是可做的
但是都学到现在了,一般我是不怎么打大整数的,而且现在可能也打不出来
有个结论,尽可能得分333,最后若是4/24/24/2,就不分了
这个结论不需要证吧,考场上肯定是找规律找出来的
可以套在FFTFFTFFT上,不就是多项式只有个常数项罢了,照样可以算
xxx个333求乘积,就对FFTFFTFFT进行快速幂,处理进位问题即可
刷怪升级——转战玄灵大陆
C:P6300
悔改
第一次见到这种应用,记在小本本上!! ✍
想了一会儿,列了个n3n^3n3的dpdpdp转移,只能搞555分,✧*。٩(㉨)و✧。 果断放弃
然后,木头拼接起来长度相加,突然想到了生成函数
发现可以用长度为iii的木头数量作为xix^ixi的系数,两端拼在一起就是生成函数的二次幂
火速开敲,马上测样例,(⊙_⊙;)?393\ 93 9怎么肥事!!
仔细想了一下,ε=(・д・`*)ハァ…两个长度为111的木头和一个长度为888的木头拼接了,算成了两根新木头,但是按题目来说必须有一个长度111木头不用
<。)#)))≦我陷入了沉思…… 🤔
我应该取两种拼接木头个数更少的一方,但是这样子怎么写呢??
到此为止,我走向了题解: 一类取min类型的卷积到普通乘法卷积的转化
∑i+j=kmin(toti,totj)=∑d=1∑i+j=k[ai≥d][aj≥d]\sum_{i+j=k}min(tot_i,tot_j)=\sum_{d=1}\sum_{i+j=k}[a_i\ge d][a_j\ge d]∑i+j=kmin(toti,totj)=∑d=1∑i+j=k[ai≥d][aj≥d]
也就是说,枚举木头个数,每次都卷积
Fd(x)=∑[ai≥d]xiF_d(x)=\sum[a_i\ge d]x^iFd(x)=∑[ai≥d]xi
ansk=12∑d=1[xk]Fd(x)2ans_k=\frac{1}{2}\sum_{d=1}[x^k]F_d(x)^2ansk=21∑d=1[xk]Fd(x)2
将tottottot数组离散化直接暴力卷积
据题解说:本质不同的FdF_dFd是O(n)O(\sqrt{n})O(n),最坏情况为1+2+...=O(D2)1+2+...=O(D^2)1+2+...=O(D2)
复杂度O(mnlogm)O(m\sqrt{n}logm)O(mnlogm)常数不大
现在才明白为什么有的题目明明没有取模,也有题解代码用的NTTNTTNTT
因为NTTNTTNTT不用敲复数那么多玩意儿
只需要取一个非常大的模数使得答案永远不可能超过即可
妙(・。・)!
#include <cmath>
#include <cstdio>
#include <iostream>
#include <algorithm>
using namespace std;
#define maxn 400100
struct complex {double x, i;complex(){}complex( double X, double I ) { x = X, i = I;}
}A[maxn];const double pi = acos( -1.0 );complex operator + ( complex a, complex b ) {return complex( a.x + b.x, a.i + b.i );
}complex operator - ( complex a, complex b ) {return complex( a.x - b.x, a.i - b.i );
}complex operator * ( complex a, complex b ) {return complex( a.x * b.x - a.i * b.i, a.x * b.i + a.i * b.x );
}int n, m, len;
int tot[maxn], r[maxn], ans[maxn], tmp[maxn];void FFT( complex *v, int opt ) {for( int i = 0;i < len;i ++ )if( i < r[i] ) swap( v[i], v[r[i]] );for( int i = 1;i < len;i <<= 1 ) {complex omega( cos( pi / i ), opt * sin( pi / i ) );for( int j = 0;j < len;j += ( i << 1 ) ) {complex w( 1, 0 );for( int k = 0;k < i;k ++, w = w * omega ) {complex x = v[j + k], y = v[j + k + i] * w;v[j + k] = x + y;v[j + k + i] = x - y;}}}
}int main() {scanf( "%d %d", &n, &m );for( int i = 1, x;i <= n;i ++ ) {scanf( "%d", &x );tot[x] ++;}len = 1; int l = 0;while( len <= ( m << 1 ) ) len <<= 1, l ++;for( int i = 0;i < len;i ++ )r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );for( int i = 1;i <= m;i ++ ) tmp[i] = tot[i];sort( tmp + 1, tmp + m + 1 );int T = unique( tmp + 1, tmp + m + 1 ) - tmp - 1;for( int i = 1;i <= T;i ++ ) {for( int j = 0;j < len;j ++ )A[j] = complex( tot[j] >= tmp[i], 0 );FFT( A, 1 );for( int j = 0;j < len;j ++ )A[j] = A[j] * A[j];FFT( A, -1 );for( int j = 0;j < len;j ++ )ans[j] += ( tmp[i] - tmp[i - 1] ) * int( A[j].x / len + 0.5 );}int maxx = 0, ansLen;for( int i = 0;i < len;i ++ )if( ans[i] / 2 > maxx ) maxx = ans[i] / 2, ansLen = i;printf( "%d %d\n", maxx, ansLen );return 0;
}
D:P3763
[TJOI2017]DNA
曾经放在后缀数组板块的题又被拉出来处决了,不,应该是来处决我了🙅
字符集大小只有{A,T,G,C}\{A,T,G,C\}{A,T,G,C},考虑分开枚举处理
等于枚举字符的设为111,反之为000
以样例‘C’‘C’‘C’为例
S0:ATCGCCCTAS_0:ATCGCCCTAS0:ATCGCCCTA
S:CTTCAS\ \ :CTTCAS :CTTCA
S0:001011100S_0:001011100S0:001011100
S:10010S\ \ :10010S :10010
假设S0S_0S0的长度为nnn,SSS的长度为mmm
对于S0[k,k+m−1]S_0[k,k+m-1]S0[k,k+m−1]与S[0,m−1]S[0,m-1]S[0,m−1]相同字符的个数则表示∑i=0m−1S0[i+k]∗S[i]\sum_{i=0}^{m-1}S_0[i+k]*S[i]∑i=0m−1S0[i+k]∗S[i]
令T[m−i−1]=S[i]T[m-i-1]=S[i]T[m−i−1]=S[i],代替SSS,式子可化为∑i=0m−1S0[i+k]∗T[m−i−1]\sum_{i=0}^{m-1}S_0[i+k]*T[m-i-1]∑i=0m−1S0[i+k]∗T[m−i−1]
发现此时是可以FFTFFTFFT卷积的,将下标看作多项式xix^ixi
∑i+j=m+k−1S0[i]∗T[j]\sum_{i+j=m+k-1}S_0[i]*T[j]∑i+j=m+k−1S0[i]∗T[j]
对于四种字符,FFTFFTFFT做四遍,加在一起就是子串与碱基序列相同字符的个数了
最后判断相差是否能控制在333以内,就知道这个子串碱基能否被改成吃藕碱基了
卷积后的xm+k−1x^{m+k-1}xm+k−1项才代表原碱基序列中[k,k+m−1][k,k+m-1][k,k+m−1]与吃藕碱基的匹配,要转化一下
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
#define maxn 400100
struct complex {double x, i;complex(){}complex( double X, double I ) {x = X, i = I;}
}A[maxn], B[maxn];const double pi = acos( -1.0 );complex operator + ( complex a, complex b ) {return complex( a.x + b.x, a.i + b.i );
}complex operator - ( complex a, complex b ) {return complex( a.x - b.x, a.i - b.i );
}complex operator * ( complex a, complex b ) {return complex( a.x * b.x - a.i * b.i, a.x * b.i + a.i * b.x );
}int len;
int r[maxn];void FFT( complex *v, int opt ) {for( int i = 0;i < len;i ++ )if( i < r[i] ) swap( v[i], v[r[i]] );for( int i = 1;i < len;i <<= 1 ) {complex omega( cos( pi / i ), opt * sin( pi / i ) );for( int j = 0;j < len;j += ( i << 1 ) ) {complex w( 1, 0 );for( int k = 0;k < i;k ++, w = w * omega ) {complex x = v[j + k], y = v[j + k + i] * w;v[j + k] = x + y;v[j + k + i] = x - y;}}}
}int T;
char S0[maxn], S[maxn], ch[4] = { 'A', 'G', 'T', 'C' };
int tot[maxn];int main() {scanf( "%d", &T );while( T -- ) {scanf( "%s %s", S0, S );memset( tot, 0, sizeof( tot ) );int n = strlen( S0 ), m = strlen( S );len = 1; int l = 0;while( len < ( n + m ) ) len <<= 1, l ++;for( int i = 0;i < len;i ++ )r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );for( int j = 0;j < 4;j ++ ) {for( int i = 0;i < n;i ++ )A[i] = complex( S0[i] == ch[j], 0 );for( int i = 0;i < m;i ++ )B[m - i - 1] = complex( S[i] == ch[j], 0 );FFT( A, 1 );FFT( B, 1 );for( int i = 0;i < len;i ++ )A[i] = A[i] * B[i];FFT( A, -1 );for( int i = m - 1;i <= m + n - 1;i ++ ) tot[i - m + 1] += int( A[i].x / len + 0.5 );for( int i = 0;i < len;i ++ )A[i] = B[i] = complex( 0, 0 );}int ans = 0;for( int i = 0;i < n - m + 1;i ++ ) //tot[i]:S0[i,i+m-1]和S[0,m-1]相同字符的个数 if( m - tot[i] <= 3 ) ans ++;printf( "%d\n", ans );}return 0;
}
E:P3321
[SDOI2015]序列统计
虽然以前做过一次,但是肯定是迷迷糊糊的,现在重新再做一遍,代码就找以前的博客吧哈哈哈哈哈,ԾㅂԾ,
相乘目前来说是超越已有工具的,考虑相乘转化为相加→\rightarrow→ 指数相加→\rightarrow→生成函数
但是好像没有哪个数学定理告诉我相乘转化为指数相加后取模答案不变的
卡在这里了(ˉ▽ˉ;)…
再考虑设dpi,jdp_{i,j}dpi,j表示选了iii个数相乘的结果取模mmm等于jjj的方案数,则有很简单的转移
f[i+1][j×kmodm]=∑k∈Sf[i][j]f[i+1][j\times k\mod m]=\sum_{k∈S}f[i][j]f[i+1][j×kmodm]=∑k∈Sf[i][j]
好像也转换不了形式使其长得像FFTFFTFFT啊凸(艹皿艹 )
尽管做过,但是还是没想出来啊/(ㄒoㄒ)/,去看一下之前写的题解
b( ̄▽ ̄)d 题目没说清楚mmm是个质数!!尽管说了我也不太想得到
指数的相加可以用mmm的原根来做底数
dpdpdp转移也稍有不同
f[i+i][j]=∑j1×j2%m=jf[i][j1]×f[i][j2]f[i+i][j]=\sum_{j_1\times j_2\% m=j}f[i][j_1]\times f[i][j_2]f[i+i][j]=∑j1×j2%m=jf[i][j1]×f[i][j2]
但是大致想了75%75\%75%左右出来,已经相比于之前完全0%0\%0%好多了
F:P5641
【CSGRound2】开拓者的卓识、
(>ρ<")推半天推的是个什么玩意儿柿子噢
这个柿子从上往下推似乎推不出什么玩意儿,推了有一会儿确实不会处理
考虑反过来从下往上推
令ansr=sumk,1,rans_r=sum_{k,1,r}ansr=sumk,1,r,考虑每个aia_iai的贡献,假设aia_iai出现了cic_ici次
显然答案表示为ansr=∑i=1nai⋅cians_r=\sum_{i=1}^na_i·c_iansr=∑i=1nai⋅ci
cic_ici还有一种含义理解:位置iii在[1,r][1,r][1,r]中被k−1k-1k−1重区间覆盖的方案数
(减一是因为最大的区间已经是被确定为[1,r][1,r][1,r]
等价于选择k−1k-1k−1个区间[l1,r1],[l2,r2]...[lk−1,rk−1][l_1,r_1],[l_2,r_2]...[l_{k-1},r_{k-1}][l1,r1],[l2,r2]...[lk−1,rk−1]
使得1≤li,ri≤r,[li,ri]⊆[li+1,ri+1],i∈[li,ri]1\le l_i,r_i\le r,[l_i,r_i]\subseteq[l_{i+1},r_{i+1}],i∈[l_i,r_i]1≤li,ri≤r,[li,ri]⊆[li+1,ri+1],i∈[li,ri]的方案数
或者…༼ ∗ •̀ (oo) •́ ∗ ༽可以枚举理解
sum1,l,rsum_{1,l,r}sum1,l,r:就是ala_lal到ara_rar的和
sum2,l,rsum_{2,l,r}sum2,l,r:可以看作选择一个区间[l1,r1]⊆[l,r][l_1,r_1]\subseteq [l,r][l1,r1]⊆[l,r],然后把al1a_{l_1}al1到ar1a_{r_1}ar1的和再加到答案里面
sum3,l,rsum_{3,l,r}sum3,l,r:选择一个区间[l2,r2]⊆[l1,r1][l_2,r_2]\subseteq [l_1,r_1][l2,r2]⊆[l1,r1],再对其做一遍sum2,l2,r2sum_{2,l_2,r_2}sum2,l2,r2,相当于再选一个区间[l2,r2]⊆[l1,r1][l_2,r_2]\subseteq [l_1,r_1][l2,r2]⊆[l1,r1],将al2a_{l_2}al2到ar2a_{r_2}ar2的和再累计到答案中
………….
sumk,l,rsum_{k,l,r}sumk,l,r:就是选择k−1k-1k−1个区间
继续等价于lil_ili在[1,i][1,i][1,i]中选,rir_iri在[i,r][i,r][i,r]中选,一共在[1,i][1,i][1,i]中选k−1k-1k−1个(可重位置,在[i,r][i,r][i,r]中选k−1k-1k−1个(可重位置,的方案数
令TnmT_{n}^mTnm表示在nnn个盒子中选择mmm个盒子的方案数
则cic_ici方案数即为Tik−1×Tr−i+1k−1T_{i}^{k-1}\times T_{r-i+1}^{k-1}Tik−1×Tr−i+1k−1
TnmT_n^mTnm也是非常好求的,就是数学中经典的隔板法
用m−1m-1m−1个板将盒子分成mmm份,每份合并看作是最后的一个超级盒子,那么一共就选择了mmm个超级盒子,符合TTT的含义
所以就可以用组合数来计算,此时就是不重位置的挑选了,Tnm=Cn+m−1mT_n^m=C_{n+m-1}^mTnm=Cn+m−1m
那么ci=C(i)+(k−1)−1k−1×C(r−i+1)+(k−1)−1k−1=Ci+k−2k−1×Cr−i+k−1k−1c_i=C_{(i)+(k-1)-1}^{k-1}\times C_{(r-i+1)+(k-1)-1}^{k-1}=C_{i+k-2}^{k-1}\times C_{r-i+k-1}^{k-1}ci=C(i)+(k−1)−1k−1×C(r−i+1)+(k−1)−1k−1=Ci+k−2k−1×Cr−i+k−1k−1
ansr=∑i=1naiCi+k−2k−1×Cr−i+k−1k−1ans_r=\sum_{i=1}^na_iC_{i+k-2}^{k-1}\times C_{r-i+k-1}^{k-1}ansr=∑i=1naiCi+k−2k−1×Cr−i+k−1k−1
令gi=Ci+k−1k−1,fi=Ci+k−2k−1g_i=C_{i+k-1}^{k-1},f_i=C_{i+k-2}^{k-1}gi=Ci+k−1k−1,fi=Ci+k−2k−1
则ansr=∑i=1ngr−i∗fians_r=\sum_{i=1}^ng_{r-i}*f_iansr=∑i=1ngr−i∗fi
就可以卷积了,有模数选择NTTNTTNTT
由于kkk非常的大,预处理组合数显然不可做,那就使用组合数的递推式即可
巧妙发现fi=ai×gi−1f_i=a_i\times g_{i-1}fi=ai×gi−1,fif_ifi也可递推获得
#include <cstdio>
#include <iostream>
using namespace std;
#define mod 998244353
#define int long long
#define maxn 400100
int n, k, len;
int a[maxn], r[maxn], f[maxn], g[maxn];int qkpow( int x, int y ) {int ans = 1;while( y ) {if( y & 1 ) ans = ans * x % mod;x = x * x % mod;y >>= 1;}return ans;
}void NTT( int *v, int opt ) {for( int i = 0;i < len;i ++ )if( i < r[i] ) swap( v[i], v[r[i]] );for( int i = 1;i < len;i <<= 1 ) {int omega = qkpow( opt == 1 ? 3 : mod / 3 + 1, ( mod - 1 ) / ( i << 1 ) );for( int j = 0;j < len;j += ( i << 1 ) )for( int k = 0, w = 1;k < i;k ++, w = w * omega % mod ) {int x = v[j + k], y = v[j + k + i] * w % mod;v[j + k] = ( x + y ) % mod;v[j + k + i] = ( x - y + mod ) % mod;}}if( opt == -1 ) {int inv = qkpow( len, mod - 2 );for( int i = 0;i < len;i ++ )v[i] = v[i] * inv % mod;}
}signed main() {scanf( "%lld %lld", &n, &k );for( int i = 1;i <= n;i ++ )scanf( "%lld", &a[i] );g[0] = 1, g[1] = k % mod, f[1] = a[1];for( int i = 2;i <= n;i ++ ) {g[i] = g[i - 1] * ( i + k - 1 ) % mod * qkpow( i, mod - 2 ) % mod;f[i] = a[i] * g[i - 1] % mod;}len = 1; int l = 0;while( len <= ( n << 1 ) ) len <<= 1, l ++;for( int i = 0;i < len;i ++ )r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );NTT( f, 1 );NTT( g, 1 );for( int i = 0;i < len;i ++ )f[i] = f[i] * g[i] % mod;NTT( f, -1 );for( int i = 1;i <= n;i ++ )printf( "%lld ", f[i] );return 0;
}
G:P4986
逃离
题意大概是这个亚子。走的范围是一个环,红线代表C(x)∗tC(x)*tC(x)∗t,绿线代表A(x)∗tA(x)*tA(x)∗t,蓝线代表B(x)∗tB(x)*tB(x)∗t
非常小学鸡化容易操作,将蓝绿线各自平移成一条线段,然后与圆的半径构成直角三角形,要从一个点一起出去,则斜边就是$ hdxrie 的走法,直角边就是的走法,直角边就是的走法,直角边就是 Althen $的走法
根据直角三角形法则,有(A(x)∗t)2+(B(x)∗t)2=(C(x)∗t)2(A(x)*t)^2+(B(x)*t)^2=(C(x)*t)^2(A(x)∗t)2+(B(x)∗t)2=(C(x)∗t)2
化简得A2(x)+B2(x)=C2(x)A^2(x)+B^2(x)=C^2(x)A2(x)+B2(x)=C2(x)
定义一个新函数f(x)=C2(x)−A2(x)−B2(x)f(x)=C^2(x)-A^2(x)-B^2(x)f(x)=C2(x)−A2(x)−B2(x)
最后的即是求函数f(x)f(x)f(x)的零点,当然必须满足x∈[L,R]x∈[L,R]x∈[L,R]
而求一个函数的根近似解,套用一阶牛顿迭代
x=x0−f(x0)f′(x0)x=x_0-\frac{f(x_0)}{f'(x_0)}x=x0−f′(x0)f(x0)
随便找一个点x0x_0x0,然后做这个点的切线,斜率就是对点x0x_0x0求导(f’(x0)f’(x_0)f’(x0))
从这个切线的根(与xxx轴的交点)出发,做一根垂线,和曲线相交于xxx点
不停重复上述过程,直到逼近原函数的零点附近
稍微解释一下牛顿迭代的公式
对于一条直线解析式y=kx+by=kx+by=kx+b,那么关于x0x_0x0的切线即为f(x0)=f′(x0)x0+bf(x_0)=f'(x_0)x_0+bf(x0)=f′(x0)x0+b
根据初中知识知道,切线斜率还可以写成y−y0x−x0\frac{y-y_0}{x-x_0}x−x0y−y0,那么用切线的零点和x0x_0x0重写切线解析式f′(x0)=0−f(x0)x−x0⇒x=x0−f(x0)f′(x0)f'(x_0)=\frac{0-f(x_0)}{x-x_0}\Rightarrow x=x_0-\frac{f(x_0)}{f'(x_0)}f′(x0)=x−x00−f(x0)⇒x=x0−f′(x0)f(x0)
最后就是设置一个牛顿迭代次数tititi,暴算去逼近零点(如果有的话
当然这过程中要和L,RL,RL,R比大小保证是答案是在范围内的
#include <cmath>
#include <cstdio>
#include <iostream>
using namespace std;
#define maxn 400100
#define eps 1e-10
struct complex {double x, i;complex(){}complex( double X, double I ) {x = X, i = I;}
}A[maxn];const double pi = acos( -1.0 );complex operator + ( complex a, complex b ) {return complex( a.x + b.x, a.i + b.i );
}complex operator - ( complex a, complex b ) {return complex( a.x - b.x, a.i - b.i );
}complex operator * ( complex a, complex b ) {return complex( a.x * b.x - a.i * b.i, a.x * b.i + a.i * b.x );
}int len;
int r[maxn];void FFT( complex *v, int opt ) {for( int i = 0;i < len;i ++ )if( i < r[i] ) swap( v[i], v[r[i]] );for( int i = 1;i < len;i <<= 1 ) {complex omega( cos( pi / i ), opt * sin( pi / i ) );for( int j = 0;j < len;j += ( i << 1 ) ) {complex w( 1, 0 );for( int k = 0;k < i;k ++, w = w * omega ) {complex x = v[j + k], y = v[j + k + i] * w;v[j + k] = x + y, v[j + k + i] = x - y;}}}
}int La, Lb, Lc, n;
double L, R;
int a[maxn], b[maxn], c[maxn], f[maxn], g[maxn];void work( int &n, int *ans ) {len = 1; int l = 0;while( len <= ( n << 1 ) ) len <<= 1, l ++;for( int i = 0;i < len;i ++ )r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );for( int i = 0;i < len;i ++ )A[i] = complex( ans[i], 0 );FFT( A, 1 );for( int i = 0;i < len;i ++ )A[i] = A[i] * A[i];FFT( A, -1 );for( int i = 0;i < len;i ++ )ans[i] = int( A[i].x / len + 0.5 );n = len;
}double calc_f( double x ) {double ans = 0, mi = 1;for( int i = 0;i <= n;i ++ )ans += f[i] * mi, mi *= x;return ans;
}double calc_g( double x ) {double ans = 0, mi = 1;for( int i = 0;i < n;i ++ )ans += g[i] * mi, mi *= x;return ans;
}int main() {scanf( "%d %d %d %lf %lf", &La, &Lb, &Lc, &L, &R );for( int i = 0;i <= La;i ++ )scanf( "%d", &a[i] );for( int i = 0;i <= Lb;i ++ )scanf( "%d", &b[i] );for( int i = 0;i <= Lc;i ++ )scanf( "%d", &c[i] );work( La, a );work( Lb, b );work( Lc, c );n = max( La, max( Lb, Lc ) );for( int i = 0;i <= n;i ++ )f[i] = c[i] - a[i] - b[i];for( int i = 0;i < n;i ++ )g[i] = f[i + 1] * ( i + 1 );int ti = 50;double x = ( L + R ) / 2.0;while( ti ) {double tmp = calc_f( x );if( fabs( tmp ) < eps ) break;x = x - tmp / calc_g( x );x = max( x, L ), x = min( x, R );ti --;}if( ti ) printf( "%.8f\n", x );else printf( "Inconsistent!" );return 0;
}
H:P4721——获得知识药剂一瓶——分治FFTFFTFFT
【模板】分治 FFT
对于蒟蒻来说还是算一个新知识,嘚学Y(>_<、)Y
分治FFT是一种基于CDQ分治的算法,主要用于在O(nlog2n)O(nlog^2n)O(nlog2n)的时间复杂度计算已知g1g_1g1到gn−1g_{n-1}gn−1,f0=1f_0=1f0=1,求fi=∑j=1ifi−j⋅gjf_i=\sum_{j=1}^if_{i-j}·g_jfi=∑j=1ifi−j⋅gj
考虑分治l,r,mid=⌊l+r2⌋l,r,mid=\lfloor\frac{l+r}{2}\rfloorl,r,mid=⌊2l+r⌋,计算左半段[l,mid][l,mid][l,mid]对于右半段[mid+1,r][mid+1,r][mid+1,r]的贡献
设x∈[mid+1,r]x∈[mid+1,r]x∈[mid+1,r],那么对fxf_xfx的贡献转移wxw_xwx即为wx=∑i=lmidfigx−iw_x=\sum_{i=l}^{mid}f_ig_{x-i}wx=∑i=lmidfigx−i
不妨将推导范围扩大一点,先保证fi=0,i∈[mid+1,r]f_i=0,i∈[mid+1,r]fi=0,i∈[mid+1,r]
wx=∑i=lx−1figx−i=∑i=0x−l−1fi+l⋅gx−l−iw_x=\sum_{i=l}^{x-1}f_ig_{x-i}=\sum_{i=0}^{x-l-1}f_{i+l}·g_{x-l-i}wx=∑i=lx−1figx−i=∑i=0x−l−1fi+l⋅gx−l−i
令Ai=fi+l,Bi=gi+1A_i=f_{i+l},B_i=g_{i+1}Ai=fi+l,Bi=gi+1,代回上式wx=∑i=0x−l−1AiBx−l−i−1w_x=\sum_{i=0}^{x-l-1}A_iB_{x-l-i-1}wx=∑i=0x−l−1AiBx−l−i−1
发现此时的A,BA,BA,B可以进行FFTFFTFFT卷积
每次操作的长度为r−l+1r-l+1r−l+1,总长度nlognnlognnlogn,一共lognlognlogn层,所以时间复杂度大概在O(nlog2n)O(nlog^2n)O(nlog2n)
#include <cstdio>
#include <iostream>
using namespace std;
#define mod 998244353
#define int long long
#define maxn 400100
int r[maxn], f[maxn], g[maxn], A[maxn], B[maxn];int qkpow( int x, int y ) {int ans = 1;while( y ) {if( y & 1 ) ans = ans * x % mod;x = x * x % mod;y >>= 1;}return ans;
}void NTT( int *v, int opt, int len ) {for( int i = 0;i < len;i ++ )if( i < r[i] ) swap( v[i], v[r[i]] );for( int i = 1;i < len;i <<= 1 ) {int omega = qkpow( opt == 1 ? 3 : mod / 3 + 1, ( mod - 1 ) / ( i << 1 ) );for( int j = 0;j < len;j += ( i << 1 ) )for( int k = 0, w = 1;k < i;k ++, w = w * omega % mod ) {int x = v[j + k], y = v[j + k + i] * w % mod;v[j + k] = ( x + y ) % mod;v[j + k + i] = ( x - y + mod ) % mod;}}if( opt == -1 ) {int inv = qkpow( len, mod - 2 );for( int i = 0;i < len;i ++ )v[i] = v[i] * inv % mod;}
}void mul( int *A, int *B, int n, int l ) {for( int i = 0;i < n;i ++ ) r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );NTT( A, 1, n );NTT( B, 1, n );for( int i = 0;i < n;i ++ ) A[i] = A[i] * B[i] % mod;NTT( A, -1, n );
}void CDQ( int L, int R ) {if( L == R ) return;int mid = ( L + R ) >> 1;CDQ( L, mid );int len = 1, l = 0;while( len <= R - L + 1 ) len <<= 1, l ++;for( int i = L;i <= mid;i ++ ) A[i - L] = f[i];for( int i = 1;i <= R - L;i ++ ) B[i - 1] = g[i];mul( A, B, len, l );for( int i = mid + 1;i <= R;i ++ )f[i] = ( f[i] + A[i - L - 1] ) % mod;for( int i = 0;i < len;i ++ ) A[i] = B[i] = 0;CDQ( mid + 1, R );
}signed main() {int n;scanf( "%lld", &n );for( int i = 1;i < n;i ++ )scanf( "%lld", &g[i] );f[0] = 1;CDQ( 0, n - 1 );for( int i = 0;i < n;i ++ )printf( "%lld ", f[i] );return 0;
}
I:P3338
[ZJOI2014]力
Fj=∑i=1j−1qi×qj(i−j)2−∑i=j+1nqi×qj(i−j)2,Ej=Fjqj⇒Ej=∑i=1j−1qi(i−j)2−∑i=j+1nqi(i−j)2F_j=\sum_{i=1}^{j-1}\frac{q_i\times q_j}{(i-j)^2}-\sum_{i=j+1}^n\frac{q_i\times q_j}{(i-j)^2},E_j=\frac{F_j}{q_j}\Rightarrow E_j=\sum_{i=1}^{j-1}\frac{q_i}{(i-j)^2}-\sum_{i=j+1}^n\frac{q_i}{(i-j)^2}Fj=i=1∑j−1(i−j)2qi×qj−i=j+1∑n(i−j)2qi×qj,Ej=qjFj⇒Ej=i=1∑j−1(i−j)2qi−i=j+1∑n(i−j)2qi
把i=ji=ji=j加进去,其实是不影响的,代码又不会真的去除,我们只一心想要卷积
∑i=1jqi(i−j)2−∑i=jnqi(i−j)2\sum_{i=1}^{j}\frac{q_i}{(i-j)^2}-\sum_{i=j}^n\frac{q_i}{(i-j)^2}∑i=1j(i−j)2qi−∑i=jn(i−j)2qi
设fi=qi,gi=1i2f_i=q_i,g_i=\frac{1}{i^2}fi=qi,gi=i21
∑i=1jfigj−i−∑i=jnfigi−j\sum_{i=1}^jf_ig_{j-i}-\sum_{i=j}^nf_ig_{i-j}∑i=1jfigj−i−∑i=jnfigi−j
令f0=g0=0f_0=g_0=0f0=g0=0
∑i=0jfi∗gj−i−∑i=jnfigi−j\sum_{i=0}^jf_i*g_{j-i}-\sum_{i=j}^nf_ig_{i-j}∑i=0jfi∗gj−i−∑i=jnfigi−j
左边已经是一个卷积的形式了,接下来只考虑右边即可
一般卷积形式都会把求和符号的下标换成从000开始
并且用两个数组的下标去凑成一个求和上限的恒定值
∑i=jnfigi−j=∑i=0n−jfi+jgi\sum_{i=j}^nf_ig_{i-j}=\sum_{i=0}^{n-j}f_{i+j}g_{i}∑i=jnfigi−j=∑i=0n−jfi+jgi
设hi=fn−ih_i=f_{n-i}hi=fn−i
∑i=0n−jhn−i−j∗gi\sum_{i=0}^{n-j}h_{n-i-j}*g_i∑i=0n−jhn−i−j∗gi
两边都变成了完美的卷积形式(๑′ᴗ‵๑)
∑i=0jf[i]∗g[j−i]−∑i=0n−jh[n−i−j]∗g[i]\sum_{i=0}^jf[i]*g[j-i]-\sum_{i=0}^{n-j}h[n-i-j]*g[i]∑i=0jf[i]∗g[j−i]−∑i=0n−jh[n−i−j]∗g[i]
设多项式F(x)=∑i=0nfi,G(x)=∑i=0ngi,H(x)=∑i=0nhiF(x)=\sum_{i=0}^nf_i,G(x)=\sum_{i=0}^ng_i,H(x)=\sum_{i=0}^nh_iF(x)=∑i=0nfi,G(x)=∑i=0ngi,H(x)=∑i=0nhi
令L(x)=F(x)∗G(x),R(x)=G(x)∗H(x)L(x)=F(x)*G(x),R(x)=G(x)*H(x)L(x)=F(x)∗G(x),R(x)=G(x)∗H(x)
则ansi=li−rn−ians_i=l_i-r_{n-i}ansi=li−rn−i,li,rn−il_i,r_{n-i}li,rn−i分别表示多项式对应项的系数
#include <cmath>
#include <cstdio>
#include <iostream>
using namespace std;
#define maxn 400100
struct complex {double x, i;complex(){}complex( double X, double I ) {x = X, i = I;}
}F[maxn], G[maxn], H[maxn];const double pi = acos( -1.0 );complex operator + ( complex a, complex b ) {return complex( a.x + b.x, a.i + b.i );
}complex operator - ( complex a, complex b ) {return complex( a.x - b.x, a.i - b.i );
}complex operator * ( complex a, complex b ) {return complex( a.x * b.x - a.i * b.i, a.x * b.i + a.i * b.x );
}int len;
int r[maxn];void FFT( complex *v, int opt ) {for( int i = 0;i < len;i ++ )if( i < r[i] ) swap( v[i], v[r[i]] );for( int i = 1;i < len;i <<= 1 ) {complex omega( cos( pi / i ), opt * sin( pi / i ) );for( int j = 0;j < len;j += ( i << 1 ) ) {complex w( 1, 0 );for( int k = 0;k < i;k ++, w = w * omega ) {complex x = v[j + k], y = v[j + k + i] * w;v[j + k] = x + y;v[j + k + i] = x - y;}}}
}int n;
int main() {scanf( "%d", &n );for( int i = 1;i <= n;i ++ ) {scanf( "%lf", &F[i].x ), F[i].i = 0;H[n - i] = F[i];G[i] = complex( 1.0 / i / i, 0 );}len = 1; int l = 0;while( len <= ( n << 1 ) ) len <<= 1, l ++;for( int i = 0;i < len;i ++ )r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );FFT( F, 1 );FFT( G, 1 );FFT( H, 1 );for( int i = 0;i < len;i ++ ) F[i] = F[i] * G[i], H[i] = H[i] * G[i];FFT( F, -1 );FFT( H, -1 );for( int i = 1;i <= n;i ++ )printf( "%.3f\n", F[i].x / len - H[n - i].x / len );return 0;
}
J:P4173
残缺的字符串
读完题,迅速反应类似的本系列的D题DNA
考虑分成262626个字母独立考虑(包括∗*∗,所在位置一直为111即可
将当前枚举的字母所在位置赋为111,其余为000
对于B[k,k+m−1]B[k,k+m-1]B[k,k+m−1]和A[0,m−1]A[0,m-1]A[0,m−1]的相同的字符长度个数表示为∑i=0m−1B[k+i]A[i]\sum_{i=0}^{m-1}B[k+i]A[i]∑i=0m−1B[k+i]A[i]
令Tm−i−1=AiT_{m-i-1}=A_{i}Tm−i−1=Ai
∑i=0m−1Bk+i∗Tm−i−1\sum_{i=0}^{m-1}B_{k+i}*T_{m-i-1}∑i=0m−1Bk+i∗Tm−i−1
同样的卷积
∑i+j=m−k−1Bi∗Tj\sum_{i+j=m-k-1}B_i*T_{j}∑i+j=m−k−1Bi∗Tj
最后toti=[m−i−1]B∗Ttot_i=[m-i-1]B*Ttoti=[m−i−1]B∗T,将272727个相同字符个数加起来,判断是否长度达到mmm即可
麻溜敲完,很不幸——303030,剩下全是TLETLETLE,╮(╯-╰)╭好吧…
你以为是再次中场休息??还是太年轻了骚年;又是一个新的好算法呢
——luoguluoguluogu第一篇题解真的很良心٩(๑`н´๑)۶
普通的单模式串匹配
长度为mmm的模式串AAA,长度为nnn的文本串BBB,求所有位置,满足BBB串从第xxx个字符开始的连续mmm个字符与AAA完全相同
Step1:定义匹配函数
定义字符串的匹配函数为C(x,y)=A(x)−B(y)C(x,y)=A(x)-B(y)C(x,y)=A(x)−B(y)
在这里形式化的定义“匹配”,若C(x,y)=0C(x,y)=0C(x,y)=0,则AAA的第xxx个字符与BBB的第yyy个字符相匹配
Step2:定义完全匹配函数
定义D(x)=∑i=0m−1C(i,x−m+i+1)D(x)=\sum_{i=0}^{m-1}C(i,x-m+i+1)D(x)=∑i=0m−1C(i,x−m+i+1),若D(x)=0D(x)=0D(x)=0,则AAA与BBB以xxx结尾的连续mmm个字符完全匹配
但到这里,发现因为带了求和符号,所以只要两串所含字符一样,势必会有D(x)=0D(x)=0D(x)=0
这就与我们的初衷背道而驰了,而导致这种情况发生的原因——C(x,y)C(x,y)C(x,y)自带±±±符号
考虑如何去掉这个符号的影响,使得只要两个串的某一位置不匹配,完全匹配函数就一定不为000
非常初中数学老师教的——加绝对值,但是在解不等式的时候,老师教过绝对值也不好弄,处理成——平方
此时有,D(x)=∑i=0m−1C(i,x−m+i+1)2=∑i=0m−1(A(i)−B(x−m+i+1))2D(x)=\sum_{i=0}^{m-1}C(i,x-m+i+1)^2=\sum_{i=0}^{m-1}\bigg(A(i)-B(x-m+i+1)\bigg)^2D(x)=∑i=0m−1C(i,x−m+i+1)2=∑i=0m−1(A(i)−B(x−m+i+1))2
Step3:Step3:Step3:暴力整式子,观察结构,不谋手段の优化
先展开看看
D(x)=∑i=0m−1(A2(i)−2A(i)B(x−m+i+1)+B2(x−m+i+1))D(x)=\sum_{i=0}^{m-1}\bigg(A^2(i)-2A(i)B(x-m+i+1)+B^2(x-m+i+1)\bigg)D(x)=∑i=0m−1(A2(i)−2A(i)B(x−m+i+1)+B2(x−m+i+1))
很套路的尝试反转AAA串,定义A′(i)=A(m−i−1)A'(i)=A(m-i-1)A′(i)=A(m−i−1)
D(x)=∑i=0m−1(A′2(m−i−1)−2A′(m−i−1)B(x−m+i+1)+B2(x−m+i+1))D(x)=\sum_{i=0}^{m-1}\bigg(A'^2(m-i-1)-2A'(m-i-1)B(x-m+i+1)+B^2(x-m+i+1)\bigg)D(x)=∑i=0m−1(A′2(m−i−1)−2A′(m−i−1)B(x−m+i+1)+B2(x−m+i+1))
一般只想求和符号里面放一个东西,求和里面的加减乘除考虑各自拆开
D(x)=∑i=0m−1A′(m−i−1)2+∑i=0m−1B(x−m+i+1)2−2∑i=0m−1A′(m−i−1)B(x−m+i+1)D(x)=\sum_{i=0}^{m-1}A'(m-i-1)^2+\sum_{i=0}^{m-1}B(x-m+i+1)^2-2\sum_{i=0}^{m-1}A'(m-i-1)B(x-m+i+1)D(x)=∑i=0m−1A′(m−i−1)2+∑i=0m−1B(x−m+i+1)2−2∑i=0m−1A′(m−i−1)B(x−m+i+1)
第一项,可以O(m)O(m)O(m)预处理,与BBB无瓜,sovledsovledsovled
第二项,可以O(n)O(n)O(n)预处理前缀和,solvedsolvedsolved
第三项,看到ta,应该是满心欢喜的,๑乛◡乛๑卷积小宝贝,改写成2∑i+j=xA′(i)∗B(j)2\sum_{i+j=x}A'(i)*B(j)2∑i+j=xA′(i)∗B(j)
最后整理一下,设T=∑i=0m−1A′(i)2,pre(x)=∑i=0xB(i)2,g(x)=∑i+j=xA′(i)∗B(j)T=\sum_{i=0}^{m-1}A'(i)^2,pre(x)=\sum_{i=0}^{x}B(i)^2,g(x)=\sum_{i+j=x}A'(i)*B(j)T=∑i=0m−1A′(i)2,pre(x)=∑i=0xB(i)2,g(x)=∑i+j=xA′(i)∗B(j)
进而有,D(x)=T+pre(x)−pre(x−m)−2g(x)D(x)=T+pre(x)-pre(x-m)-2g(x)D(x)=T+pre(x)−pre(x−m)−2g(x)
FFTFFTFFT优化算g(x)g(x)g(x),最后O(n)O(n)O(n)算D(x)D(x)D(x)
通配符的单模式串匹配
考虑依葫芦画瓢,一样的步骤
Step1:定义匹配函数
设通配符的数值为000,000乘任何数为000,既然想把与通配符的匹配算出来的C(x,y)C(x,y)C(x,y)也变成000,不放把数值乘进去
C(x,y)=(A(x)−B(y))2⋅A(x)⋅B(y)C(x,y)=\big(A(x)-B(y)\big)^2·A(x)·B(y)C(x,y)=(A(x)−B(y))2⋅A(x)⋅B(y)
Step2:完全匹配函数
D(x)=∑i=0m−1(A(i)−B(x−m+i+1))2⋅A(i)⋅B(x−m+i+1)D(x)=\sum_{i=0}^{m-1}\big(A(i)-B(x-m+i+1)\big)^2·A(i)·B(x-m+i+1)D(x)=∑i=0m−1(A(i)−B(x−m+i+1))2⋅A(i)⋅B(x−m+i+1)
Step3:暴力化式子
D(x)=∑i=0m−1A(i)3B(x−m+i+1)+∑i=0m−1B(x−m+i+1)3A(i)−2∑i=0m−1A(i)2B(x−m+i+1)2D(x)=\sum_{i=0}^{m-1}A(i)^3B(x-m+i+1)+\sum_{i=0}^{m-1}B(x-m+i+1)^3A(i)-2\sum_{i=0}^{m-1}A(i)^2B(x-m+i+1)^2D(x)=∑i=0m−1A(i)3B(x−m+i+1)+∑i=0m−1B(x−m+i+1)3A(i)−2∑i=0m−1A(i)2B(x−m+i+1)2
定义A′(i)=A(m−i−1)A'(i)=A(m-i-1)A′(i)=A(m−i−1)
D(x)=∑i=0m−1A′(m−i−1)3B(x−m+i+1)+∑i=0m−1B(x−m+i+1)3A′(m−i−1)−2∑i=0m−1A′(x−m−1)2B(x−m+i+1)2D(x)=\sum_{i=0}^{m-1}A'(m-i-1)^3B(x-m+i+1)+\sum_{i=0}^{m-1}B(x-m+i+1)^3A'(m-i-1)-2\sum_{i=0}^{m-1}A'(x-m-1)^2B(x-m+i+1)^2D(x)=∑i=0m−1A′(m−i−1)3B(x−m+i+1)+∑i=0m−1B(x−m+i+1)3A′(m−i−1)−2∑i=0m−1A′(x−m−1)2B(x−m+i+1)2
更妙的地方出现了!!发现所有小括号加起来都是xxx!!!卷积宝贝作法了ԅ(¯㉨¯ԅ)
D(x)=∑i+j=xA′(i)3B(j)+∑i+j=xB(j)3A′(i)−2∑i+j=xA′(i)2B(j)2D(x)=\sum_{i+j=x}A'(i)^3B(j)+\sum_{i+j=x}B(j)^3A'(i)-2\sum_{i+j=x}A'(i)^2B(j)^2D(x)=∑i+j=xA′(i)3B(j)+∑i+j=xB(j)3A′(i)−2∑i+j=xA′(i)2B(j)2
三次多项式乘法即可!!!(〜㉨)〜(▽)
常数稍稍大了点,但这不重要,吸氧能过的都不叫卡常p(# ̄▽ ̄#)o
#include <cmath>
#include <cstdio>
#include <iostream>
#include <algorithm>
using namespace std;
#define maxn 1200100
struct complex {double x, i;complex(){}complex( double X, double I ) {x = X, i = I;}
}a[maxn], b[maxn], D[maxn];const double pi = acos( -1.0 );complex operator + ( complex a, complex b ) {return complex( a.x + b.x, a.i + b.i );
}complex operator - ( complex a, complex b ) {return complex( a.x - b.x, a.i - b.i );
}complex operator * ( complex a, complex b ) {return complex( a.x * b.x - a.i * b.i, a.x * b.i + a.i * b.x );
}int len;
int r[maxn];void FFT( complex *v, int opt ) {for( int i = 0;i < len;i ++ )if( i < r[i] ) swap( v[i], v[r[i]] );for( int i = 1;i < len;i <<= 1 ) {complex omega( cos( pi / i ), opt * sin( pi / i ) );for( int j = 0;j < len;j += ( i << 1 ) ) {complex w( 1, 0 );for( int k = 0;k < i;k ++, w = w * omega ) {complex x = v[j + k], y = v[j + k + i] * w;v[j + k] = x + y;v[j + k + i] = x - y;}}}
}int m, n, cnt;
char s1[maxn], s2[maxn];
int ans[maxn], A[maxn], B[maxn];int main() {scanf( "%d %d %s %s", &m, &n, s1, s2 );len = 1; int l = 0;while( len < m + n ) len <<= 1, l ++;for( int i = 0;i < len;i ++ )r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );for( int i = 0;i < m;i ++ ) A[i] = ( s1[i] == '*' ) ? 0 : s1[i] - 'a' + 1;for( int i = 0;i < n;i ++ )B[i] = ( s2[i] == '*' ) ? 0 : s2[i] - 'a' + 1;reverse( A, A + m );for( int i = 0;i < len;i ++ ) {a[i] = complex( A[i] * A[i] * A[i], 0 );b[i] = complex( B[i], 0 );}FFT( a, 1 ), FFT( b, 1 );for( int i = 0;i < len;i ++ )D[i] = D[i] + a[i] * b[i];for( int i = 0;i < len;i ++ ) {a[i] = complex( A[i], 0 );b[i] = complex( B[i] * B[i] * B[i], 0 );}FFT( a, 1 ), FFT( b, 1 );for( int i = 0;i < len;i ++ )D[i] = D[i] + a[i] * b[i];for( int i = 0;i < len;i ++ ) {a[i] = complex( A[i] * A[i], 0 );b[i] = complex( B[i] * B[i], 0 );}FFT( a, 1 ), FFT( b, 1 );for( int i = 0;i < len;i ++ )D[i] = D[i] - a[i] * b[i] * complex( 2, 0 );FFT( D, -1 );for( int i = m - 1;i < n;i ++ )if( fabs( D[i].x / len ) < 0.5 ) ans[++ cnt] = i - m + 2;//本来D[i]算的是i-m+1 但此题要求下标从1开始 所以要多加1printf( "%d\n", cnt );for( int i = 1;i <= cnt;i ++ )printf( "%d ", ans[i] );return 0;
}
K:P5748
集合划分计数
题意:求n≤1e5n\le 1e5n≤1e5的贝尔数,bell(n)bell(n)bell(n)即将nnn个有标号的球划分若干集合的方案数
设fi,jf_{i,j}fi,j表示把iii个数划分为jjj个集合的方案数
那么就有两种转移,要么iii自成一个集合,要么加入jjj个当中的任意一个集合
fi,j=fi−1,j−1+fi−1,j×jf_{i,j}=f_{i-1,j-1}+f_{i-1,j}\times jfi,j=fi−1,j−1+fi−1,j×j
ans=∑i=1nfn,ians=\sum_{i=1}^nf_{n,i}ans=∑i=1nfn,i
O(n2)O(n^2)O(n2)简单的暴力DPDPDP转移
设fif_ifi表示把nnn个数划分为iii个集合的方案数
(╯▽╰ ) 卡住了——( ̄△ ̄;)
法一:生成函数
指数型生成函数用来解决排列问题,普通型生成函数用来解决组合问题
设一个非空集合的指数型生成函数F(x)F(x)F(x),答案的指数型生成函数为G(x)G(x)G(x)
枚举划分的集合,集合间不区分,有G(x)=∑i=0∞Fi(x)i!=eF(x)G(x)=\sum_{i=0}^∞\frac{F^i(x)}{i!}=e^{F(x)}G(x)=∑i=0∞i!Fi(x)=eF(x)
一个非空集合的指数型生成函数,显然为ex−1e^x-1ex−1(∑i=0∞xii!\sum_{i=0}^∞\frac{x^i}{i!}∑i=0∞i!xi减去一个空集f0=1f_0=1f0=1
多项式expexpexp即可,题目求的是有标号的方案数,最后乘以n!n!n!即可
法二:分治FFTFFTFFT
枚举第一个元素所在集合大小,列出状态转移方程
fi=∑j=1iCi−1,j−1×fi−j=∑j=1i(i−1)!(j−1)!(i−j)!×fi−jf_i=\sum_{j=1}^iC_{i-1,j-1}\times f_{i-j}=\sum_{j=1}^i\frac{(i-1)!}{(j-1)!(i-j)!}\times f_{i-j}fi=∑j=1iCi−1,j−1×fi−j=∑j=1i(j−1)!(i−j)!(i−1)!×fi−j
把带i,ji,ji,j的项分开
fi(i−1)!=∑j=1i1(j−1)!×fi−j(i−j)!\frac{f_i}{(i-1)!}=\sum_{j=1}^i\frac{1}{(j-1)!}\times \frac{f_{i-j}}{(i-j)!}(i−1)!fi=∑j=1i(j−1)!1×(i−j)!fi−j
令Fi=1i!,Gi=fii!F_i=\frac{1}{i!},G_i=\frac{f_i}{i!}Fi=i!1,Gi=i!fi,卷积!!
fi(i−1)!=∑j=0i−1Fj∗Gi−j−1=∑j+k=i−1Fj∗Gk\frac{f_i}{(i-1)!}=\sum_{j=0}^{i-1}F_{j}*G_{i-j-1}=\sum_{j+k=i-1}F_j*G_k(i−1)!fi=∑j=0i−1Fj∗Gi−j−1=∑j+k=i−1Fj∗Gk
FFF是已知多项式,预处理即可,GGG是一个只与fj,j<if_j,j<ifj,j<i有关的多项式,最后答案别忘记乘(i−1)!(i-1)!(i−1)!
这不是左边对右边造成贡献的好套路吗??——中场休息的分治FFTFFTFFT登场!
#include <cstdio>
#include <iostream>
using namespace std;
#define int long long
#define mod 998244353
#define maxn 400100
#define n 100000
int r[maxn], fac[maxn], inv[maxn], f[maxn], g[maxn], A[maxn], B[maxn];int qkpow( int x, int y ) {int ans = 1;while( y ) {if( y & 1 ) ans = ans * x % mod;x = x * x % mod;y >>= 1;}return ans;
}void NTT( int *v, int opt, int len ) {for( int i = 0;i < len;i ++ )if( i < r[i] ) swap( v[i], v[r[i]] );for( int i = 1;i < len;i <<= 1 ) {int omega = qkpow( opt == 1 ? 3 : mod / 3 + 1, ( mod - 1 ) / ( i << 1 ) );for( int j = 0;j < len;j += ( i << 1 ) )for( int k = 0, w = 1;k < i;k ++, w = w * omega % mod ) {int x = v[j + k], y = v[j + k + i] * w % mod;v[j + k] = ( x + y ) % mod;v[j + k + i] = ( x - y + mod ) % mod; }}if( opt == -1 ) {int Inv = qkpow( len, mod - 2 );for( int i = 0;i < len;i ++ )v[i] = v[i] * Inv % mod;}
}void init() {fac[0] = inv[0] = 1;for( int i = 1;i <= n;i ++ ) fac[i] = fac[i - 1] * i % mod;inv[n] = qkpow( fac[n], mod - 2 );for( int i = n - 1;i;i -- ) inv[i] = inv[i + 1] * ( i + 1 ) % mod;
}void CDQ( int L, int R ) {if( L == R ) { if( L ) f[L] = f[L] * fac[L - 1] % mod; return; }int mid = ( L + R ) >> 1;CDQ( L, mid );int len = 1, l = 0;while( len <= R - L + 1 ) len <<= 1, l ++;for( int i = 0;i < len;i ++ ) r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );for( int i = L;i <= mid;i ++ ) A[i - L] = f[i] * inv[i] % mod;for( int i = 0;i < len;i ++ ) B[i] = inv[i];NTT( A, 1, len );NTT( B, 1, len );for( int i = 0;i < len;i ++ ) A[i] = A[i] * B[i] % mod;NTT( A, -1, len );for( int i = mid + 1;i <= R;i ++ ) f[i] = ( f[i] + A[i - 1 - L] ) % mod;for( int i = 0;i < len;i ++ ) A[i] = B[i] = 0;CDQ( mid + 1, R );
}signed main() {init();f[0] = 1;CDQ( 0, n );int T;scanf( "%lld", &T );while( T -- ) {int x;scanf( "%lld", &x );printf( "%lld\n", f[x] );}return 0;
}
L:P3702
[SDOI2017]序列计数
至少有一个为质数(;¬_¬)俺讨厌至少。。。
转化为满足要求的序列总数减去一个质数都没有的满足要求的序列数
一个质数都没有的序列数,不妨把mmm以内的所有质数全都丢掉
设fi,jf_{i,j}fi,j:前iii个数的总和取模ppp为jjj的方案数(一个质数都没有,枚举第iii个选的数
fi,j=∑(j1+j2)%p=jfi−1,j2f_{i,j}=\sum_{(j_1+j_2)\%p=j}f_{i-1,j_2}fi,j=∑(j1+j2)%p=jfi−1,j2,j1j_1j1为合数
下面带取模,好熟悉,EEE题的序列统计⊙﹏⊙!,修改转移方程
fi<<1,j=∑(j1+j2)%p=jfi,j1∗fi,j2f_{i<<1,j}=\sum_{(j_1+j_2)\%p=j}f_{i,j_1}*f_{i,j_2}fi<<1,j=∑(j1+j2)%p=jfi,j1∗fi,j2
突然想起,ppp好像没保证是质数。。。淦((٩(//̀Д/́/)۶))没有原根,指数上的取模就不能进行了
不!٩(๑`н´๑)۶我是sb,这比序列统计还要简单,下面的条件是加又不是乘,干嘛甩成指数做,多此一举!!
已经是个卷积形式了。。。难道三模数NTTNTTNTT??害怕码力.<{=....(嘎嘎嘎~)
有矩阵快速幂+容斥的,容斥实在对于现在的我要求太高了
暴力卷积快速幂就可以跑过去(•̀⌄•́乁,这算不算搞出这道题了p(# ̄▽ ̄#)o
#include <cstdio>
#define int long long
#define mod 20170408
#define maxn 20000005
#define maxp 400
int n, m, p, cnt;
int f[maxp], g[maxp], c[maxp], F[maxp], G[maxp], prime[maxn];
bool vis[maxn];void mul( int *a, int *b ) {for( int i = 0;i < p;i ++ )for( int j = 0;j < p;j ++ )c[i + j] = ( c[i + j] + a[i] * b[j] % mod ) % mod;for( int i = 0;i < p;i ++ )a[i] = ( c[i] + c[i + p] ) % mod, c[i] = c[i + p] = 0;
}signed main() {scanf( "%lld %lld %lld", &n, &m, &p );f[1] = g[1] = F[0] = G[0] = 1;for( int i = 2;i <= m;i ++ ) {f[i % p] ++;if( ! vis[i] ) prime[++ cnt] = i;else g[i % p] ++;for( int j = 1;j <= cnt && i * prime[j] <= m;j ++ ) {vis[i * prime[j]] = 1;if( i % prime[j] == 0 ) break;}}while( n ) {if( n & 1 ) mul( F, f ), mul( G, g );mul( f, f ), mul( g, g );n >>= 1;}printf( "%lld\n", ( F[0] - G[0] + mod ) % mod );return 0;
}
这题的难度放错位置了,就当中场休息好了
M:P5667
拉格朗日插值2
套拉格朗日差值公式fm+x=∑i=0nf(i)∏i≠jm+x−ji−jf_{m+x}=\sum_{i=0}^nf(i)\prod_{i≠j}\frac{m+x-j}{i-j}fm+x=∑i=0nf(i)∏i=ji−jm+x−j
∏i≠jm+x−j→(m+x)(m+x−1)+...(m+x−n)/(m+x−i)=(m+x)!(m+x−n−1)!(m+x−i)\prod_{i≠j}m+x-j\rightarrow (m+x)(m+x-1)+...(m+x-n)/(m+x-i)=\frac{(m+x)!}{(m+x-n-1)!(m+x-i)}∏i=jm+x−j→(m+x)(m+x−1)+...(m+x−n)/(m+x−i)=(m+x−n−1)!(m+x−i)(m+x)!
对于枚举的iii,有[0,i)[0,i)[0,i)小于iii,相减为正,而对于(i,n](i,n](i,n]的共n−in-in−i个数而言,分母为负
所以乘积的符号便由n−in-in−i的奇偶决定,(−1)n−i(-1)^{n-i}(−1)n−i
∏j<ii−j=(i−0)(i−1)(i−2)...(2)(1)=i!\prod_{j<i}i-j=(i-0)(i-1)(i-2)...(2)(1)=i!∏j<ii−j=(i−0)(i−1)(i−2)...(2)(1)=i!
∏i<ji−j=(n−i)(n−i+1)...(2)(1)=(n−i)!\prod_{i<j}i-j=(n-i)(n-i+1)...(2)(1)=(n-i)!∏i<ji−j=(n−i)(n−i+1)...(2)(1)=(n−i)!(不考虑正负,上面已经判断了
⇒fm+x=∑i=0n(f(i)×(m+x)!(−1)n−ii!(n−i)!(m+x−n−1)!(m+x−i))\Rightarrow f_{m+x}=\sum_{i=0}^n\bigg(f(i)\times \frac{(m+x)!}{(-1)^{n-i}i!(n-i)!(m+x-n-1)!(m+x-i)}\bigg)⇒fm+x=∑i=0n(f(i)×(−1)n−ii!(n−i)!(m+x−n−1)!(m+x−i)(m+x)!)
把只与xxx有关的式子提出来,(−1)n−i(-1)^{n-i}(−1)n−i只带来符号改变,不妨放在分子上(分母太肥了(;¬_¬)
⇒fm+x=(m+x)!(m+x−n−1)!×∑i=0n((−1)n−if(i)i!(n−i)!∗1m+x−i)\Rightarrow f_{m+x}=\frac{(m+x)!}{(m+x-n-1)!}\times \sum_{i=0}^n\bigg(\frac{(-1)^{n-i}f(i)}{i!(n-i)!}*\frac{1}{m+x-i}\bigg)⇒fm+x=(m+x−n−1)!(m+x)!×∑i=0n(i!(n−i)!(−1)n−if(i)∗m+x−i1)
发现是个卷积形式,但是iii的取值不是0≤i≤x0\le i\le x0≤i≤x,而达到了nnn,与一般的有所区别(・。・)
但是没关系啊ԅ(¯㉨¯ԅ),我们卷积最喜欢的不就是重新定义新函数改变原多项式顺序吗?只要最后能很快速地反映所求答案位置
老子把整个序列硬往后挪nnn位,第x+nx+nx+n位计算的答案存的就是原来fm+xf_{m+x}fm+x
令Ai=(−1)n−if(i)i!(n−i)!,Bi=1m−n+iA_i=\frac{(-1)^{n-i}f(i)}{i!(n-i)!},B_i=\frac{1}{m-n+i}Ai=i!(n−i)!(−1)n−if(i),Bi=m−n+i1
fx+m=(x+m)!(m+x−n−1)!∑i=0nAi∗Bn+x−i=(x+m)!(m+x−n−1)!∑j+k=n+xAj∗Bkf_{x+m}=\frac{(x+m)!}{(m+x-n-1)!}\sum_{i=0}^nA_i*B_{n+x-i}=\frac{(x+m)!}{(m+x-n-1)!}\sum_{j+k=n+x}A_j*B_kfx+m=(m+x−n−1)!(x+m)!∑i=0nAi∗Bn+x−i=(m+x−n−1)!(x+m)!∑j+k=n+xAj∗Bk
令i=x+ni=x+ni=x+n(x∈[0,n]⇒i∈[n,2n]x∈[0,n]\Rightarrow i∈[n,2n]x∈[0,n]⇒i∈[n,2n]
=(i+m−n)!(i+m−n−n−1)!∑j+k=iAj∗Bk=\frac{(i+m-n)!}{(i+m-n-n-1)!}\sum_{j+k=i}A_j*B_k=(i+m−n−n−1)!(i+m−n)!∑j+k=iAj∗Bk
预处理阶乘faci=i!fac_i=i!faci=i!,iii逆元inviinv_iinvi,i!i!i!阶乘逆元invfaciinvfac_iinvfaci,差位的阶乘MS_faci=∏j=0i−1(m−n+j)MS\_fac_i=\prod_{j=0}^{i-1}(m-n+j)MS_faci=∏j=0i−1(m−n+j),m−n+im-n+im−n+i差位逆元MS_invi+1MS\_inv_{i+1}MS_invi+1,MS_faciMS\_fac_iMS_faci差位阶乘逆元MS_facinviMS\_facinv_iMS_facinvi,则有1(i+m−n)!=MS_faci+1,(i+m−n−1−n)!=MS_facinvi−n\frac{1}{(i+m-n)!}=MS\_fac_{i+1},(i+m-n-1-n)!=MS\_facinv_{i-n}(i+m−n)!1=MS_faci+1,(i+m−n−1−n)!=MS_facinvi−n
最后bb这么多,是为了解释代码的一些地方,避免不必要的时间投入看代码,FFT就是这种玩下标的东西,经常容易玩懵,找准哪个位置是自己原来求解数组的对应位置是关键
#include <cstdio>
#include <iostream>
using namespace std;
#define int long long
#define mod 998244353
#define maxn 1000000
int n, m, len;
int f[maxn], r[maxn], A[maxn], B[maxn];
int inv[maxn], facinv[maxn], fac[maxn];
int MS_fac[maxn], MS_inv[maxn], MS_facinv[maxn];int qkpow( int x, int y ) {int ans = 1;while( y ) {if( y & 1 ) ans = ans * x % mod;x = x * x % mod;y >>= 1;}return ans;
}void NTT( int *v, int opt ) {for( int i = 0;i < len;i ++ )if( i < r[i] ) swap( v[i], v[r[i]] );for( int i = 1;i < len;i <<= 1 ) {int omega = qkpow( opt == 1 ? 3 : mod / 3 + 1, ( mod - 1 ) / ( i << 1 ) );for( int j = 0;j < len;j += ( i << 1 ) )for( int k = 0, w = 1;k < i;k ++, w = w * omega % mod ) {int x = v[j + k], y = v[j + k + i] * w % mod;v[j + k] = ( x + y ) % mod;v[j + k + i] = ( x - y + mod ) % mod;}}if( opt == -1 ) {int Inv = qkpow( len, mod - 2 );for( int i = 0;i < len;i ++ )v[i] = v[i] * Inv % mod;}
}void init( int n, int m ) {fac[0] = MS_fac[0] = MS_inv[0] = inv[0] = 1;for( int i = 1;i <= ( n << 1 | 1 );i ++ ) {fac[i] = fac[i - 1] * i % mod;MS_fac[i] = MS_fac[i - 1] * ( m - n + i - 1 ) % mod;}facinv[n << 1 | 1] = qkpow( fac[n << 1 | 1], mod - 2 );MS_facinv[n << 1 | 1] = qkpow( MS_fac[n << 1 | 1], mod - 2 );for( int i = ( n << 1 );~ i;i -- ) {facinv[i] = facinv[i + 1] * ( i + 1 ) % mod;MS_facinv[i] = MS_facinv[i + 1] * ( m - n + i ) % mod;}for( int i = ( n << 1 | 1 );i;i -- ) {inv[i] = facinv[i] * fac[i - 1] % mod;MS_inv[i] = MS_facinv[i] * MS_fac[i - 1] % mod;}
}signed main() {scanf( "%lld %lld", &n, &m );init( n, m );for( int i = 0;i <= n;i ++ )scanf( "%lld", &f[i] );len = 1; int l = 0;while( len <= n * 3 ) len <<= 1, l ++;for( int i = 0;i < len;i ++ )r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );for( int i = 0;i <= n;i ++ ) {A[i] = f[i] * facinv[i] % mod * facinv[n - i] % mod;if( ( n - i ) & 1 ) A[i] = mod - A[i];}for( int i = 0;i <= ( n << 1 );i ++ )B[i] = MS_inv[i + 1];NTT( A, 1 ), NTT( B, 1 );for( int i = 0;i < len;i ++ )A[i] = A[i] * B[i] % mod;NTT( A, -1 );for( int i = n;i <= ( n << 1 );i ++ )printf( "%lld ", MS_fac[i + 1] * MS_facinv[i - n] % mod * A[i] % mod );return 0;
}
N:P4245——拾得知识药丸一颗——三模数NTTNTTNTT
【模板】任意模数多项式乘法
别称:三模数NTTNTTNTT 嘿嘿嘿——才发现这种三模数NTTNTTNTT就是MTTMTTMTT,那没事了V●ᴥ●V
当然这题有很多好的解法,但是鄙人认为还是应该掌握基础的不要脑子的sb基础做法
{x≡x1(modM1)x≡x2(modM2)x≡x3(modM3)\begin{cases}x\equiv x_1\ \ (mod\ \ M_1)\\x\equiv x_2\ \ (mod\ \ M_2)\\x\equiv x_3\ \ (mod\ \ M_3)\end{cases}⎩⎪⎨⎪⎧x≡x1 (mod M1)x≡x2 (mod M2)x≡x3 (mod M3)
x=x1+k1M1=x2+k2M2⇒x1+k1M1≡x2(modM2)⇒k1=x2−x1M1(modM2)x=x_1+k_1M_1=x_2+k_2M_2\Rightarrow x_1+k_1M_1\equiv x_2\pmod {M_2}\Rightarrow k_1=\frac{x_2-x_1}{M_1}\pmod {M_2}x=x1+k1M1=x2+k2M2⇒x1+k1M1≡x2(modM2)⇒k1=M1x2−x1(modM2)
x≡x1+k1M1(modM1M2)x\equiv x_1+k_1M_1\pmod {M_1M_2}x≡x1+k1M1(modM1M2),令t=x1+k1M1t=x_1+k_1M_1t=x1+k1M1
t+k′M1M2=x3+k3M3⇒t+k′M1M2≡x3(modM3)⇒k′=x3−tM1M2(modM3)t+k'M_1M_2=x_3+k_3M_3\Rightarrow t+k'M_1M_2\equiv x_3\pmod{M_3}\Rightarrow k'=\frac{x_3-t}{M_1M_2}\pmod{M_3}t+k′M1M2=x3+k3M3⇒t+k′M1M2≡x3(modM3)⇒k′=M1M2x3−t(modM3)
∵x<M1M2M3⇒x=k′M1M2+t∵x<M_1M_2M_3\Rightarrow x=k'M_1M_2+t∵x<M1M2M3⇒x=k′M1M2+t
中国剩余定理及扩展剩余定理证明
#include <cstdio>
#include <iostream>
using namespace std;
#define int long long
#define maxn 400100
int len, n, m, p;
int r[maxn], ans[maxn], B[maxn], F[maxn], G[maxn];
int Mod[5] = { 0, 998244353, 469762049, 1004535809 };int qkpow( int x, int y, int mod ) {x %= mod;int ret = 1;while( y ) {if( y & 1 ) ret = ret * x % mod;x = x * x % mod;y >>= 1;}return ret;
}struct poly {int A[maxn], mod;void NTT( int *v, int opt ) {for( int i = 0;i < len;i ++ )if( i < r[i] ) swap( v[i], v[r[i]] );for( int i = 1;i < len;i <<= 1 ) {int omega = qkpow( opt == 1 ? 3 : mod / 3 + 1, ( mod - 1 ) / ( i << 1 ), mod );for( int j = 0;j < len;j += ( i << 1 ) )for( int k = 0, w = 1;k < i;k ++, w = w * omega % mod ) {int x = v[j + k], y = v[j + k + i] * w % mod;v[j + k] = ( x + y ) % mod;v[j + k + i] = ( x - y + mod ) % mod;}}if( opt == -1 ) {int inv = qkpow( len, mod - 2, mod );for( int i = 0;i < len;i ++ )v[i] = v[i] * inv % mod;}}
}ntt[5];int qkmul( int x, int y, int mod ) {x %= mod;int ret = 0;while( y ) {if( y & 1 ) ret = ( ret + x ) % mod;x = ( x + x ) % mod;y >>= 1;}return ret;
}void CRT() {int mod_12 = Mod[1] * Mod[2];int inv1 = qkpow( Mod[1], Mod[2] - 2, Mod[2] );int inv2 = qkpow( Mod[2], Mod[1] - 2, Mod[1] );int inv_12 = qkpow( mod_12, Mod[3] - 2, Mod[3] );for( int i = 0;i < len;i ++ ) {int x1 = ntt[1].A[i];int x2 = ntt[2].A[i];int x3 = ntt[3].A[i];int t = ( qkmul( x1 * Mod[2], inv2, mod_12 ) + qkmul( x2 * Mod[1], inv1, mod_12 ) ) % mod_12;int k = ( x3 - t % Mod[3] + Mod[3] ) % Mod[3] * inv_12 % Mod[3];ans[i] = ( k % p * ( mod_12 % p ) % p + t % p ) % p;}
}signed main() {scanf( "%lld %lld %lld", &n, &m, &p );for( int i = 0;i <= n;i ++ ) scanf( "%lld", &F[i] );for( int i = 0;i <= m;i ++ ) scanf( "%lld", &G[i] );len = 1;int l = 0;while( len <= n + m ) len <<= 1, l ++;for( int i = 0;i < len;i ++ )r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );for( int k = 1;k <= 3;k ++ ) {ntt[k].mod = Mod[k];for( int i = 0;i <= n;i ++ ) ntt[k].A[i] = F[i];for( int i = 0;i <= m;i ++ ) B[i] = G[i];for( int i = m + 1;i < len;i ++ ) B[i] = 0;ntt[k].NTT( ntt[k].A, 1 );ntt[k].NTT( B, 1 );for( int i = 0;i < len;i ++ ) ntt[k].A[i] = ntt[k].A[i] * B[i] % Mod[k];ntt[k].NTT( ntt[k].A, -1 );}CRT();for( int i = 0;i <= n + m;i ++ )printf( "%lld ", ans[i] );return 0;
}
O:P5488
差分与前缀和
一.前缀和
si=∑j=0iajs_i=\sum_{j=0}^ia_jsi=∑j=0iaj,考虑前缀和的式子与卷积式子si=∑j+k=iaj∗bks_i=\sum_{j+k=i}a_j*b_ksi=∑j+k=iaj∗bk
发现前缀和相当于是乘了一个系数全是111的多项式∑i=0∞1xi\sum_{i=0}^∞1\ x_i∑i=0∞1 xi
由等比数列公式可得该多项式的闭形式11−x\frac{1}{1-x}1−x1
kkk阶前缀和,就相当于卷了kkk次,卷积有结合律,所以最后相当于原来的数组卷积1(1−x)k\frac{1}{(1-x)^k}(1−x)k1
这个式子的第nnn项系数为Cn+k−1n=Cn+k−1k−1C_{n+k-1}^{n}=C_{n+k-1}^{k-1}Cn+k−1n=Cn+k−1k−1,递推即可
至于为什么是这个系数,可以打表验证,其实发现前缀和是跟杨辉三角有着某种联系的,意会也可
二.差分
差分相当于是乘了一个系数全是−1-1−1的多项式∑i=0∞−1xi\sum_{i=0}^∞ -1\ x^i∑i=0∞−1 xi,闭形式为1−x1-x1−x
最后相当于卷积(1−x)k=∑i=0∞(−1)iCkixi(1-x)^k=\sum_{i=0}^∞(-1)^iC_k^ix^i(1−x)k=∑i=0∞(−1)iCkixi
递推即可Cki=Cki−1×k−i+1iC_k^i=C_k^{i-1}\times \frac{k-i+1}{i}Cki=Cki−1×ik−i+1
#include <cstdio>
#include <iostream>
using namespace std;
#define mod 1004535809
#define int long long
#define maxn 400100
int len;
int r[maxn], inv[maxn], A[maxn], B[maxn];int read() {int x = 0;char s = getchar();while( s < '0' || s > '9' ) s = getchar();while( '0' <= s && s <= '9' ) {x = ( ( x << 1 ) + ( x << 3 ) + ( s - 48 ) ) % mod;s = getchar();}return x;
}int qkpow( int x, int y ) {int ans = 1;while( y ) {if( y & 1 ) ans = ans * x % mod;x = x * x % mod;y >>= 1;}return ans;
}void NTT( int *v, int opt ) {for( int i = 0;i < len;i ++ )if( i < r[i] ) swap( v[i], v[r[i]] );for( int i = 1;i < len;i <<= 1 ) {int omega = qkpow( opt == 1 ? 3 : mod / 3 + 1, ( mod - 1 ) / ( i << 1 ) );for( int j = 0;j < len;j += ( i << 1 ) ) for( int k = 0, w = 1;k < i;k ++, w = w * omega % mod ) {int x = v[j + k], y = v[j + k + i] * w % mod;v[j + k] = ( x + y ) % mod;v[j + k + i] = ( x - y + mod ) % mod;}}if( opt == -1 ) {int Inv = qkpow( len, mod - 2 );for( int i = 0;i < len;i ++ )v[i] = v[i] * Inv % mod;}
}signed main() {int n, type, k;scanf( "%lld", &n );k = read();scanf( "%lld", &type );for( int i = 0;i < n;i ++ ) scanf( "%lld", &A[i] );inv[1] = 1;for( int i = 2;i < n;i ++ ) inv[i] = ( mod - mod / i ) * inv[mod % i] % mod;B[0] = 1;if( ! type ) {for( int i = 1;i < n;i ++ )B[i] = B[i - 1] * ( i + k - 1 ) % mod * inv[i] % mod;}else {for( int i = 1;i < n;i ++ )B[i] = ( B[i - 1] * ( k - i + 1 + mod ) % mod ) % mod * inv[i] % mod;for( int i = 1;i < n;i ++ )if( i & 1 ) B[i] = mod - B[i];}len = 1;int l = 0;while( len < ( n << 1 ) ) len <<= 1, l ++;for( int i = 0;i < len;i ++ )r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );NTT( A, 1 ), NTT( B, 1 );for( int i = 0;i < len;i ++ ) A[i] = A[i] * B[i] % mod;NTT( A, -1 );for( int i = 0;i < n;i ++ ) printf( "%lld ", A[i] );return 0;
}
P:P4199
万径人踪灭
问:不连续的关于某根对称轴对称的回文子序列个数
非常难做,正面硬刚刚不动 (-ロ-)(-ロ-)
考虑转化问题:关于某根对称轴对称的回文子序列个数-回文子串个数
而对于一个序列求其回文子串个数,直接运用manachermanachermanacher即可求解
现在考虑怎么求关于某根对称轴对称的回文子序列个数
对称轴有两种情况,假设左右对称两个位置字符相等个数为kkk(不算对称轴
①以某个位置为轴
答案为2k+1−12^{k+1}-12k+1−1,那些位置以及对称轴选与不选再减去全都不选的情况
②以某两个位置的间隙为轴
答案为2k−12^k-12k−1
考虑每个位置,fif_ifi:以iii为对称轴的两个位置字符相等的个数
fi=∑j=0i[sj,si×2−j]f_i=\sum_{j=0}^i[s_j,s_{i\times 2-j}]fi=∑j=0i[sj,si×2−j]
非常可以,这很卷积,套用DDD题DNADNADNA的套路,分字符自己卷自己,再把个数加起来
注意:其实卷积后求出的fif_ifi表示以i2\frac{i}{2}2i为对称轴(除不尽就是以位置间隙为轴的意思
位置字符相等个数×2\times 2×2
所以真正的个数应该为fi2\frac{f_i}{2}2fi
但是如果i2\frac{i}{2}2i除得尽,代表一个位置的话,卷积中的gi2∗gi2g_\frac{i}{2}*g_{\frac{i}{2}}g2i∗g2i便只算了一次,
不妨用⌊fi2⌋+[2∣i]\lfloor\frac{f_i}{2}\rfloor+[2|i]⌊2fi⌋+[2∣i]来特别加上即可,统计答案就是∑i=02n−22⌊fi2⌋+[2∣i]−1\sum_{i=0}^{2n-2}2^{\lfloor\frac{f_i}{2}\rfloor+[2|i]}-1∑i=02n−22⌊2fi⌋+[2∣i]−1
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
#define maxn 400100
#define int long long
#define MOD 1000000007
#define mod 998244353
int len;
int rev[maxn], A[maxn], B[maxn], f[maxn];int qkpow( int x, int y, int Mod ) {int ans = 1;while( y ) {if( y & 1 ) ans = ans * x % Mod;x = x * x % Mod;y >>= 1;}return ans;
}void NTT( int *v, int opt ) {for( int i = 0;i < len;i ++ )if( i < rev[i] ) swap( v[i], v[rev[i]] );for( int i = 1;i < len;i <<= 1 ) {int omega = qkpow( opt == 1 ? 3 : mod / 3 + 1, ( mod - 1 ) / ( i << 1 ), mod );for( int j = 0;j < len;j += ( i << 1 ) )for( int k = 0, w = 1;k < i;k ++, w = w * omega % mod ) {int x = v[j + k], y = v[j + k + i] * w % mod;v[j + k] = ( x + y ) % mod;v[j + k + i] = ( x - y + mod ) % mod;}}if( opt == -1 ) {int inv = qkpow( len, mod - 2, mod );for( int i = 0;i < len;i ++ )v[i] = v[i] * inv % mod;}
}int r[maxn], p[maxn];
char s[maxn];int manacher() {int n = ( strlen( s ) << 1 | 1 ), idx = 0;for( int i = 0;i < n;i ++ )p[i] = ( i & 1 ) ? s[idx ++] : '#';int R = -1, C = -1;for( int i = 0;i < n;i ++ ) {r[i] = ( R > i ) ? min( R - i, r[C * 2 - i] ) : 1;while( i + r[i] < n && ~ ( i - r[i] ) )if( p[i + r[i]] == p[i - r[i]] ) r[i] ++;else break;if( i + r[i] > R ) R = i + r[i], C = i;}int ans = 0;for( int i = 0;i < n;i ++ )ans += ( r[i] >> 1 );return ans % MOD;
}signed main() {scanf( "%s", s );int n = strlen( s );len = 1;int l = 0;while( len < ( n << 1 ) ) len <<= 1, l ++;for( int i = 0;i < len;i ++ )rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );for( int i = 0;i < n;i ++ )if( s[i] == 'a' ) A[i] = 1;for( int i = 0;i < n;i ++ )if( s[i] == 'b' ) B[i] = 1;NTT( A, 1 ), NTT( B, 1 );for( int i = 0;i < len;i ++ )f[i] = ( A[i] * A[i] % mod + B[i] * B[i] % mod ) % mod;NTT( f, -1 );int ans = 0;for( int i = 0, k = 1;i <= n * 2 - 2;i ++, k ^= 1 )ans = ( ans + qkpow( 2, ( f[i] >> 1 ) + k, MOD ) - 1 + MOD ) % MOD;printf( "%lld\n", ( ans - manacher() + MOD ) % MOD );return 0;
}
Q:P4061
[HEOI2016/TJOI2016]求和
fn=∑i=0n∑j=0iS(i,j)×2j×j!f_n=\sum_{i=0}^n\sum_{j=0}^iS(i,j)\times 2^j\times j!fn=∑i=0n∑j=0iS(i,j)×2j×j!
∵S(i,j)=0[i<j]∵ S(i,j)=0[i<j]∵S(i,j)=0[i<j]
∴fn=∑i=0n∑j=0nS(i,j)×2j×j!∴f_n=\sum_{i=0}^n\sum_{j=0}^nS(i,j)\times 2^j\times j!∴fn=∑i=0n∑j=0nS(i,j)×2j×j!
S(n,m)=1m!∑i=0m(−1)iCmi(m−i)n=∑i=0n(−1)ii!(m−i)n(m−i)!S(n,m)=\frac{1}{m!}\sum_{i=0}^m(-1)^iC_m^i(m-i)^n=\sum_{i=0}^n\frac{(-1)^i}{i!}\frac{(m-i)^n}{(m-i)!}S(n,m)=m!1∑i=0m(−1)iCmi(m−i)n=∑i=0ni!(−1)i(m−i)!(m−i)n
fn=∑j=0n2j×j!∑i=0n∑k=0j(−1)kk!⋅(j−k)i(j−k)!=∑j=0n2j×j!∑k=0j(−1)kk!⋅∑i=0n(j−k)i(j−k)!f_n=\sum_{j=0}^n2^j\times j!\sum_{i=0}^n\sum_{k=0}^j\frac{(-1)^k}{k!}·\frac{(j-k)^i}{(j-k)!}=\sum_{j=0}^n2^j\times j!\sum_{k=0}^j\frac{(-1)^k}{k!}·\frac{\sum_{i=0}^n(j-k)^i}{(j-k)!}fn=∑j=0n2j×j!∑i=0n∑k=0jk!(−1)k⋅(j−k)!(j−k)i=∑j=0n2j×j!∑k=0jk!(−1)k⋅(j−k)!∑i=0n(j−k)i
(换一下变量名更好看V●ᴥ●V
fn=∑i=0n2i×i!∑j=0i(−1)jj!⋅∑k=0n(i−j)k(i−j)!f_n=\sum_{i=0}^n2^i\times i!\sum_{j=0}^i\frac{(-1)^j}{j!}·\frac{\sum_{k=0}^n(i-j)^k}{(i-j)!}fn=∑i=0n2i×i!∑j=0ij!(−1)j⋅(i−j)!∑k=0n(i−j)k
令gi=(−1)ii!,hi=∑k=0niki!=in+1−1(i−1)i!g_i=\frac{(-1)^i}{i!},h_i=\frac{\sum_{k=0}^ni^k}{i!}=\frac{i^{n+1}-1}{(i-1)i!}gi=i!(−1)i,hi=i!∑k=0nik=(i−1)i!in+1−1(等比数列公式不可能还有人呢不会吧嘿嘿嘿qwq
fn=∑i=0n2i×i!∑j=0ngj∗hi−jf_n=\sum_{i=0}^n2^i\times i!\sum_{j=0}^ng_j*h_{i-j}fn=∑i=0n2i×i!∑j=0ngj∗hi−j(特别的,g0=0,g1=n+1g_0=0,g_1=n+1g0=0,g1=n+1
NTTNTTNTT冲鸭!!
#include <cstdio>
#include <iostream>
using namespace std;
#define int long long
#define mod 998244353
#define maxn 400100
int n, m, len;
int r[maxn], f[maxn], g[maxn];
int fac[maxn], inv[maxn], facinv[maxn];int qkpow( int x, int y ) {int ans = 1;while( y ) {if( y & 1 ) ans = ans * x % mod;x = x * x % mod;y >>= 1;}return ans;
}void NTT( int *v, int opt ) {for( int i = 0;i < len;i ++ )if( i < r[i] ) swap( v[i], v[r[i]] );for( int i = 1;i < len;i <<= 1 ) {int omega = qkpow( opt == 1 ? 3 : mod / 3 + 1, ( mod - 1 ) / ( i << 1 ) );for( int j = 0;j < len;j += ( i << 1 ) )for( int k = 0, w = 1;k < i;k ++, w = w * omega % mod ) {int x = v[j + k], y = v[j + k + i] * w % mod;v[j + k] = ( x + y ) % mod;v[j + k + i] = ( x - y + mod ) % mod;}} if( opt == -1 ) {int Inv = qkpow( len, mod - 2 );for( int i = 0;i < len;i ++ )v[i] = v[i] * Inv % mod;}
}void init( int n ) {fac[0] = inv[0] = inv[1] = facinv[0] = 1;for( int i = 1;i <= n;i ++ )fac[i] = fac[i - 1] * i % mod;facinv[n] = qkpow( fac[n], mod - 2 );for( int i = n - 1;i;i -- )facinv[i] = facinv[i + 1] * ( i + 1 ) % mod;for( int i = 2;i <= n;i ++ )inv[i] = ( mod - mod / i ) * inv[mod % i] % mod;
}signed main() {scanf( "%lld", &n );init( n );for( int i = 0;i <= n;i ++ ) {g[i] = facinv[i];if( i & 1 ) g[i] = mod - g[i];}f[0] = 1, f[1] = n + 1;for( int i = 2;i <= n;i ++ )f[i] = ( qkpow( i, n + 1 ) - 1 ) * inv[i - 1] % mod * facinv[i] % mod;len = 1; int l = 0;while( len <= ( n << 1 ) ) len <<= 1, l ++;for( int i = 0;i < len;i ++ )r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );NTT( g, 1 ), NTT( f, 1 ); for( int i = 0;i < len;i ++ )f[i] = f[i] * g[i] % mod;NTT( f, -1 );int ans = 0;for( int i = 0;i <= n;i ++ )ans = ( ans + qkpow( 2, i ) * fac[i] % mod * f[i] % mod ) % mod;printf( "%lld\n", ans );return 0;
}
缘遇逆天秘籍——金牌打野
搭配生成函数修炼事半功倍哦♪(´ε`)
R:P4841
[集训队作业2013]城市规划
求nnn个点的简单有标号无向连通图数目
法一:分治FFTFFTFFT
设fif_ifi表示iii个点的有标号无向连通图数,考虑列状态转移方程
呃—— (-ロ-) 好像转移不动
再设gig_igi表示iii个点的图的数量,有i×(i−1)2\frac{i\times (i-1)}{2}2i×(i−1)条边,每条可选可不选,很容易就得出gi=2i×(i−1)2g_i=2^\frac{i\times (i-1)}{2}gi=22i×(i−1)
gn=∑i=1nCn−1i−1fign−ig_n=\sum_{i=1}^nC_{n-1}^{i-1}f_ig_{n-i}gn=∑i=1nCn−1i−1fign−i
枚举111号点所在连通图的大小,联通的方案数即为fif_ifi,有标号所以从剩下的点中选i−1i-1i−1个,乘Cn−1i−1C_{n-1}^{i-1}Cn−1i−1,其余点分成什么样的 图都可以,方案数为gn−ig_{n-i}gn−i
gn=∑i=1nCn−1i−1fign−i=Cn−1n−1fng0+∑i=1n−1Cn−1i−1fign−i=fn+∑i=1n−1(n−1)!(i−1)!(n−i)!fign−ig_n=\sum_{i=1}^nC_{n-1}^{i-1}f_ig_{n-i}=C_{n-1}^{n-1}f_ng_0+\sum_{i=1}^{n-1}C_{n-1}^{i-1}f_ig_{n-i}=f_n+\sum_{i=1}^{n-1}\frac{(n-1)!}{(i-1)!(n-i)!}f_ig_{n-i}gn=∑i=1nCn−1i−1fign−i=Cn−1n−1fng0+∑i=1n−1Cn−1i−1fign−i=fn+∑i=1n−1(i−1)!(n−i)!(n−1)!fign−i
⇒fn(n−1)!=gn(n−1)!−∑i=1n−1fi(i−1)!gn−i(n−i)!\Rightarrow \frac{f_n}{(n-1)!}=\frac{g_n}{(n-1)!}-\sum_{i=1}^{n-1}\frac{f_i}{(i-1)!}\frac{g_{n-i}}{(n-i)!}⇒(n−1)!fn=(n−1)!gn−∑i=1n−1(i−1)!fi(n−i)!gn−i
令Fi=fi(i−1)!,Gi=gii!F_i=\frac{f_i}{(i-1)!},G_i=\frac{g_i}{i!}Fi=(i−1)!fi,Gi=i!gi,则Fn=Gn×n−∑i=1n−1Fi∗Gn−iF_n=G_n\times n-\sum_{i=1}^{n-1}F_i*G_{n-i}Fn=Gn×n−∑i=1n−1Fi∗Gn−i
分治FFTFFTFFT卷积宝贝又出现了mua!(*╯3╰)
#include <cstdio>
#include <iostream>
using namespace std;
#define mod 1004535809
#define int long long
#define maxn 520100
int len;
int r[maxn], fac[maxn], inv[maxn], f[maxn], g[maxn], F[maxn], G[maxn];int qkpow( int x, int y ) {int ans = 1;while( y ) {if( y & 1 ) ans = ans * x % mod;x = x * x % mod;y >>= 1;}return ans;
}void NTT( int *v, int opt, int len ) {for( int i = 0;i < len;i ++ )if( i < r[i] ) swap( v[i], v[r[i]] );for( int i = 1;i < len;i <<= 1 ) {int omega = qkpow( opt == 1 ? 3 : mod / 3 + 1, ( mod - 1 ) / ( i << 1 ) );for( int j = 0;j < len;j += ( i << 1 ) )for( int k = 0, w = 1;k < i;k ++, w = w * omega % mod ) {int x = v[j + k], y = v[j + k + i] * w % mod;v[j + k] = ( x + y ) % mod;v[j + k + i] = ( x - y + mod ) % mod;}}if( opt == -1 ) {int Inv = qkpow( len, mod - 2 );for( int i = 0;i < len;i ++ )v[i] = v[i] * Inv % mod;}
}void solve( int L, int R ) {if( L == R ) return;int mid = ( L + R ) >> 1;solve( L, mid );int len = 1, l = 0;while( len <= R - L + 1 ) len <<= 1, l ++;for( int i = 0;i < len;i ++ )r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );for( int i = L;i <= mid;i ++ ) F[i - L] = f[i];for( int i = 1;i <= R - L;i ++ ) G[i - 1] = g[i];NTT( F, 1, len ), NTT( G, 1, len );for( int i = 0;i < len;i ++ ) F[i] = F[i] * G[i] % mod;NTT( F, -1, len );for( int i = mid + 1;i <= R;i ++ ) f[i] = ( f[i] - F[i - L - 1] + mod ) % mod;for( int i = 0;i < len;i ++ ) F[i] = G[i] = 0;solve( mid + 1, R );
}void init( int n ) {fac[0] = inv[0] = 1;for( int i = 1;i <= n;i ++ ) fac[i] = fac[i - 1] * i % mod;inv[n] = qkpow( fac[n], mod - 2 );for( int i = n - 1;i;i -- ) inv[i] = inv[i + 1] * ( i + 1 ) % mod;for( int i = 1;i <= n;i ++ )g[i] = qkpow( 2, i * ( i - 1 ) / 2 ) * inv[i] % mod, f[i] = g[i] * i % mod;
}signed main() {int n;scanf( "%lld", &n );init( maxn - 1 );solve( 0, n );printf( "%lld\n", f[n] * fac[n - 1] % mod );return 0;
}
法二:生成函数+多项式求逆
gn=∑i=1nCn−1i−1fign−ig_n=\sum_{i=1}^nC_{n-1}^{i-1}f_ig_{n-i}gn=∑i=1nCn−1i−1fign−i
换一种推式子的方式
gn=∑i=1n(n−1)!(i−1)!(n−i)!fign−i⇒gn(n−1)!=∑i=1nfi(i−1)!gn−i(n−i)!g_n=\sum_{i=1}^n\frac{(n-1)!}{(i-1)!(n-i)!}f_ig_{n-i}\Rightarrow \frac{g_n}{(n-1)!}=\sum_{i=1}^{n}\frac{f_i}{(i-1)!}\frac{g_{n-i}}{(n-i)!}gn=∑i=1n(i−1)!(n−i)!(n−1)!fign−i⇒(n−1)!gn=∑i=1n(i−1)!fi(n−i)!gn−i
令G′,F,GG',F,GG′,F,G分别为它们的生成函数,则G′(x)=∑i=0∞gi+1xii!,F(x)=∑i=0∞fi+1xii!,G(x)=∑i=0∞gixii!G'(x)=\sum_{i=0}^∞g_{i+1}\frac{x^i}{i!},F(x)=\sum_{i=0}^∞f_{i+1}\frac{x^i}{i!},G(x)=\sum_{i=0}^∞g_i\frac{x^i}{i!}G′(x)=∑i=0∞gi+1i!xi,F(x)=∑i=0∞fi+1i!xi,G(x)=∑i=0∞gii!xi
G′=F∗G⇒F=G′GG'=F*G\Rightarrow F=\frac{G'}{G}G′=F∗G⇒F=GG′
多项式求逆即可
#include <cmath>
#include <cstdio>
#include <iostream>
using namespace std;
#define int long long
#define mod 1004535809
#define maxn 600000
int n, len, l;
int r[maxn], a[maxn], b[maxn], c[maxn], d[maxn], fac[maxn];int qkpow( int x, int y ) {int ans = 1;while( y ) {if( y & 1 ) ans = ans * x % mod;x = x * x % mod;y >>= 1;}return ans;
}void NTT( int *c, int f ) {for( int i = 0;i < len;i ++ ) if( i < r[i] ) swap( c[i], c[r[i]] );for( int i = 1;i < len;i <<= 1 ) {int omega = qkpow( f == 1 ? 3 : 334845270, ( mod - 1 ) / ( i << 1 ) );for( int j = 0;j < len;j += ( i << 1 ) ) {int w = 1;for( int k = 0;k < i;k ++, w = w * omega % mod ) {int x = c[j + k], y = c[j + k + i] * w % mod;c[j + k] = ( x + y ) % mod;c[j + k + i] = ( x - y + mod ) % mod;}}}if( f == -1 ) {int inv = qkpow( len, mod - 2 );for( int i = 0;i < len;i ++ ) c[i] = c[i] * inv % mod;}
}void polyinv( int deg, int *b ) {if( deg == 1 ) { b[0] = 1; return; }polyinv( deg + 1 >> 1, b );for( len = 1, l = 0;len < ( deg << 1 ) - 1;len <<= 1, l ++ );for( int i = 0;i < len;i ++ ) r[i] = r[i >> 1] >> 1 | ( ( i & 1 ) << l - 1 );for( int i = 0;i < deg;i ++ ) c[i] = a[i];for( int i = deg;i < n;i ++ ) c[i] = 0;NTT( c, 1 ), NTT( b, 1 );for( int i = 0;i < len;i ++ ) b[i] = ( 2 - b[i] * c[i] % mod + mod ) % mod * b[i] % mod;NTT( b, -1 );for( int i = deg;i < n;i ++ ) b[i] = 0;
}signed main() {scanf( "%lld", &n ); fac[0] = 1, n ++;for( int i = 1;i <= n;i ++ ) fac[i] = fac[i - 1] * i % mod;a[0] = 1;for( int i = 1;i < n;i ++ ) {a[i] = qkpow( 2, i * ( i - 1 ) / 2 ) * qkpow( fac[i], mod - 2 ) % mod;d[i] = qkpow( 2, i * ( i - 1 ) / 2 ) * qkpow( fac[i - 1], mod - 2 ) % mod;}polyinv( n, b );for( len = 1, l = 0;len < ( n << 1 );len <<= 1, l ++ );for( int i = 0;i < len;i ++ ) r[i] = r[i >> 1] >> 1 | ( ( i & 1 ) << l - 1 );NTT( d, 1 ), NTT( b, 1 );for( int i = 0;i < len;i ++ ) d[i] = d[i] * b[i] % mod;NTT( d, -1 );printf( "%lld\n", d[n - 1] * fac[n - 2] % mod );return 0;
}
法三:生成函数+多项式lnlnln
设fif_ifi表示iii个点的无标号的无向连通图数量,gig_igi表示无标号的图数量,令F(x),G(x)F(x),G(x)F(x),G(x)为别为其的EGFEGFEGF,枚举连通块个数
G(x)=∑i=0∞F(x)ii!=eF(x)⇒F(x)=lnG(x)G(x)=\sum_{i=0}^∞\frac{F(x)^i}{i!}=e^{F(x)}\Rightarrow F(x)=lnG(x)G(x)=∑i=0∞i!F(x)i=eF(x)⇒F(x)=lnG(x)
最后[xn][x^n][xn]系数即为答案,因为题目是有标号的连通图数,所以要乘n!n!n!,同理gig_igi是无标号的数,G(x)G(x)G(x)的系数应为2i×(i−1)2i!\frac{2^\frac{i\times (i-1)}{2}}{i!}i!22i×(i−1)
S:[UVA4671]K-neighbor substrings
读完题目发现这是 D 的弱化版。同样的思路,枚举字符设为 111 其余为 000,卷积求出该字符匹配数量。
所有字符匹配数量求和,然后判断是否 ≥m−K\ge m-K≥m−K 即可。
唯一的不同就是,这个是求 本质不同的子串 数量,首先想到的就是哈希。
但是很离谱,单模双模哈希怎么样都会挂,unsigned long long
自然溢出就可以通过。🙄
#include <map>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
#define ull unsigned long long
#define maxn 800005
map < ull, bool > mp;
int len, l;
int r[maxn], ans[maxn];
ull Hash[maxn], mi[maxn];
char a[maxn], b[maxn];struct complex {double x, i;complex(){}complex( double X, double I ) { x = X, i = I; }
}f[maxn], g[maxn], h[maxn];
double pi = acos( -1.0 );
complex operator + ( complex u, complex v ) {return complex( u.x + v.x, u.i + v.i );
}
complex operator - ( complex u, complex v ) {return complex( u.x - v.x, u.i - v.i );
}
complex operator * ( complex u, complex v ) {return complex( u.x * v.x - u.i * v.i, u.i * v.x + u.x * v.i );
}void FFT( complex *c, int opt ) {for( int i = 0;i < len;i ++ ) if( i < r[i] ) swap( c[i], c[r[i]] );for( int i = 1;i < len;i <<= 1 ) {complex omega( cos( pi / i ), opt * sin( pi / i ) );for( int j = 0;j < len;j += ( i << 1 ) ) {complex w( 1, 0 );for( int k = 0;k < i;k ++, w = w * omega ) {complex x = c[j + k], y = c[j + k + i] * w;c[j + k] = x + y, c[j + k + i] = x - y;}}}
}bool check( int l, int r ) {ull num = Hash[r] - Hash[l - 1] * mi[r - l + 1];if( mp[num] ) return 0;else { mp[num] = 1; return 1; }
}int main() {mi[0] = 1;for( int i = 1;i < maxn;i ++ ) mi[i] = mi[i - 1] * 131;int K, Case = 0;while( scanf( "%d", &K ) and ~ K ) {scanf( "%s %s", a, b );int n = strlen( a ), m = strlen( b );for( len = 1, l = 0;len < ( n + m );len <<= 1, l ++ ) len <<= 1, l ++;for( int i = 0;i < len;i ++ ) r[i] = r[i >> 1] >> 1 | ( i & 1 ) << l - 1; for( char ch = 'a';ch <= 'b';ch ++ ) {for( int i = 0;i < len;i ++ ) {f[i] = complex( ch == a[i], 0 );h[i] = complex( ch == b[i], 0 );g[i] = complex( 0, 0 );}for( int i = 0;i < m;i ++ ) g[i] = h[m - i - 1];FFT( f, 1 ), FFT( g, 1 );for( int i = 0;i < len;i ++ ) f[i] = f[i] * g[i];FFT( f, -1 );for( int i = m - 1;i < len;i ++ ) ans[i - m + 1] += (int)( f[i].x / len + 0.5 );}Hash[0] = a[0] - 'a';for( int i = 1;i < n;i ++ ) Hash[i] = Hash[i - 1] * 131 + (a[i] - 'a');int ret = 0;for( int i = 0;i <= n - m;i ++ ) if( m - ans[i] <= K and check( i, i + m - 1 ) ) ret ++;printf( "Case %d: %d\n", ++ Case, ret );for( int i = 0;i < len;i ++ ) ans[i] = 0;mp.clear();}return 0;
}
T:SP8372 TSUM - Triple Sums
洛谷链接
Ai+Aj+Ak=SA_i+A_j+A_k=SAi+Aj+Ak=S 已经属于是条件反射了,将值当成下标,直接卷积就行 ansAi+Aj+Ak=Sans_{A_i+A_j+A_k=S}ansAi+Aj+Ak=S 。
但是这里有两个问题,一个是忽略了 i<j<ki<j<ki<j<k 的顺序要求,在最后答案 /6/6/6 即可。另一个是可能算重,因为 AAA 并没有互不相同。
算重就容斥好了。B2Ai,C3AiB_{2A_i},C_{3A_i}B2Ai,C3Ai。
B,CB,CB,C 同样考虑顺序问题,所以最后的 ansi=ai3−3aibi+2cians_{i}=a_i^3-3a_ib_i+2c_iansi=ai3−3aibi+2ci。再 IDFT\text{IDFT}IDFT 回去。
注意 AAA 有负的可能,直接整体平移 200002000020000 即可,这些都是小细节。
#include <cmath>
#include <cstdio>
#include <iostream>
using namespace std;
#define int long long
#define maxn 500000
int len, l;
int r[maxn];struct complex {double x, i;complex(){}complex( double X, double I ) { x = X; i = I; }
}a[maxn], b[maxn], c[maxn], ans[maxn];
double pi = acos( -1.0 );
complex operator + ( complex u, complex v ) {return complex( u.x + v.x, u.i + v.i );
}
complex operator - ( complex u, complex v ) {return complex( u.x - v.x, u.i - v.i );
}
complex operator * ( complex u, complex v ) {return complex( u.x * v.x - u.i * v.i, u.x * v.i + u.i * v.x );
}void FFT( complex *c, int opt ) {for( int i = 0;i < len;i ++ ) if( i < r[i] ) swap( c[i], c[r[i]] );for( int i = 1;i < len;i <<= 1 ) {complex omega( cos( pi / i ), opt * sin( pi / i ) );for( int j = 0;j < len;j += ( i << 1 ) ) {complex w( 1, 0 );for( int k = 0;k < i;k ++, w = w * omega ) {complex x = c[j + k], y = w * c[j + k + i];c[j + k] = x + y, c[j + k + i] = x - y;}}}
}signed main() {int n; scanf( "%lld", &n );for( int i = 1, x;i <= n;i ++ ) {scanf( "%lld", &x ); x += 20000;a[x].x ++;b[x * 2].x ++;c[x * 3].x ++;}for( len = 1, l = 0;len <= 150000;len <<= 1, l ++ );for( int i = 0;i < len;i ++ ) r[i] = r[i >> 1] >> 1 | (i & 1) << l - 1;FFT( a, 1 ), FFT( b, 1 ), FFT( c, 1 );for( int i = 0;i < len;i ++ ) ans[i] = a[i] * a[i] * a[i] - complex( 3, 0 ) * a[i] * b[i] + complex( 2, 0 ) * c[i];FFT( ans, -1 );for( int i = 0;i < len;i ++ ) {int x = (int)( ans[i].x / len + 0.5 );x /= 6;if( x > 0 ) printf( "%lld : %lld\n", i - 60000, x );}return 0;
}
U:[CodeChef - COUNTARI]Arithmetic Progressions
分块和 FFTFFTFFT。
选择一个适当的块大小 blockblockblock,然后枚举每个块。
将合法的三元组 (Ai,Aj,Ak)(A_i,A_j,A_k)(Ai,Aj,Ak) 分为两种情况考虑。
- 在块内。
枚举块内的两个点,一个做 AiA_iAi,一个做 Aj/AkA_j/A_kAj/Ak。
如果第三个是后面块的,就需要知道后面未操作块中是该值的个数。
如果第三个是前面块或者本块的,因为本块内元素的枚举有序所以不会算重,就需要知道前面已操作块和本块已枚举 AiA_iAi 的数是该值的个数。
分别用 suf,pre,cntsuf,pre,cntsuf,pre,cnt 来记录后面未操作块,前面已操作块和本块已操作数的信息。 - 在块外。
AjA_jAj 在本块内,AiA_iAi 在已操作块中,AkA_kAk 在未操作块中。
这就是很常见的卷积了。这也是为什么上一种情况需要单独记录本块已操作数量,而不是直接与 preprepre 合并。
最后每个块统计答案的时候,顺便将 cnt→precnt\rightarrow precnt→pre 中。
选择好合适的块大小,假设值最大值为 mmm,时间复杂度为 O(n2block+block⋅mlogm)O(\frac{n^2}{block}+block·m\log m)O(blockn2+block⋅mlogm)。
#include <cmath>
#include <cstdio>
#include <iostream>
using namespace std;
#define maxn 200000
int len, l, n, block, Max;
long long ans;
int r[maxn], a[maxn], pre[maxn], cnt[maxn], suf[maxn];struct complex {double x, i;complex(){}complex( double X, double I ) { x = X, i = I; }
}f[maxn], g[maxn];
double pi = acos( -1.0 );
complex operator + ( complex u, complex v ) {return complex( u.x + v.x, u.i + v.i );
}
complex operator - ( complex u, complex v ) {return complex( u.x - v.x, u.i - v.i );
}
complex operator * ( complex u, complex v ) {return complex( u.x * v.x - u.i * v.i, u.x * v.i + u.i * v.x );
}void FFT( complex *c, int opt ) {for( int i = 0;i < len;i ++ ) if( i < r[i] ) swap( c[i], c[r[i]] );for( int i = 1;i < len;i <<= 1 ) {complex omega( cos( pi / i ), opt * sin( pi / i ) );for( int j = 0;j < len;j += (i << 1) ) {complex w( 1, 0 );for( int k = 0;k < i;k ++, w = w * omega ) {complex x = c[j + k], y = c[j + k + i] * w;c[j + k] = x + y, c[j + k + i] = x - y;}}}if( opt == -1 ) for( int i = 0;i < len;i ++ ) c[i].x /= len;
}int main() {scanf( "%d", &n );block = sqrt( n );for( int i = 1;i <= n;i ++ ) {scanf( "%d", &a[i] );Max = max( Max, a[i] );suf[a[i]] ++;}for( len = 1;len <= (Max << 1);len <<= 1, l ++ );for( int i = 0;i < len;i ++ ) r[i] = r[i >> 1] >> 1 | (i & 1) << l - 1;for( int l = 1;l <= n;l += block ) {int r = min( l + block - 1, n ), k;for( int i = l;i <= r;i ++ ) suf[a[i]] --;for( int i = l;i <= r;i ++ ) {for( int j = i + 1;j <= r;j ++ ) {k = (a[i] << 1) - a[j];if( 1 <= k and k <= Max ) ans += cnt[k] + pre[k];k = (a[j] << 1) - a[i];if( 1 <= k and k <= Max ) ans += suf[k];}cnt[a[i]] ++;}for( int i = 0;i <= Max;i ++ ) {f[i] = complex( pre[i], 0 );g[i] = complex( suf[i], 0 );}for( int i = Max + 1;i < len;i ++) {f[i] = complex( 0, 0 );g[i] = complex( 0, 0 );}FFT( f, 1 ), FFT( g, 1 );for( int i = 0;i < len;i ++ ) f[i] = f[i] * g[i]; FFT( f, -1 );for( int i = l;i <= r;i ++ ) {ans += (long long)( f[a[i] << 1].x + 0.5 );pre[a[i]] ++, cnt[a[i]] --;}}printf( "%lld\n", ans );return 0;
}
V:CodeForces528D Fuzzy Search
洛谷链接
做了很多个这种 {A,T,G,C}\text{\{A,T,G,C\}}{A,T,G,C} 的字符 FFT/NTT\text{FFT/NTT}FFT/NTT 的题后,我发现 这种多个字符之间的适配问题,一般都是拆开单个单个字符做的。
考虑最基础的适配是要求完全相同,即 k=0k=0k=0 的严苛条件,我们套路地采取的是翻转一个串,gi′=gm−i−1g'_i=g_{m-i-1}gi′=gm−i−1,使得卷积下标和为定值,分字符卷积后求和,判断是否和数组长度相同。
这是解决字符适配问题的通解基础,所有题目都只是在这个的基础上进行系列的改动罢了。
对于一个字符 chchch。
设 fi=1/0f_i=1/0fi=1/0 表示 SiS_iSi 偏差 kkk 以内是否有字符 chchch 出现,即 fi=[∃j∈[i−k+1,i+k−1]Sj=ch]f_i=\big[\exist_{j\in[i-k+1,i+k-1]\ }S_j=ch\big]fi=[∃j∈[i−k+1,i+k−1] Sj=ch]。
设 gi=1/0g_i=1/0gi=1/0 表示 TiT_iTi 是否为字符 chchch,即 gi=[Ti=ch]g_i=\big[T_i=ch\big]gi=[Ti=ch]。
同样地,翻转 ggg,gi′=gm−i−1g'_{i}=g_{m-i-1}gi′=gm−i−1,然后和 fif_ifi 卷积后为 hih_ihi。
则,如果 hi≠0h_i\neq 0hi=0 也就意味着有 hih_ihi 个 Ta=chT_a=chTa=ch 在 kkk 偏差内有对应的 Sb=chS_b=chSb=ch。
注意:通过题意可以知道有可能 Ta1,Ta2T_{a_1},T_{a_2}Ta1,Ta2 对应的是同一个 bbb。
多个字符累和 hih_ihi 后,对于 hi=mh_i=mhi=m 的 iii 位置就意味着整个 TTT 都适配成功。
#include <cmath>
#include <cstdio>
#include <iostream>
using namespace std;
#define maxn 1000000
int r[maxn];
int len, l;struct complex {double x, i;complex(){}complex( double X, double I ) { x = X, i = I; }
}f[maxn], g[maxn];
double pi = acos( -1.0 );
complex operator + ( complex u, complex v ) {return complex( u.x + v.x, u.i + v.i );
}
complex operator - ( complex u, complex v ) {return complex( u.x - v.x, u.i - v.i );
}
complex operator * ( complex u, complex v ) {return complex( u.x * v.x - u.i * v.i, u.x * v.i + u.i * v.x );
}void FFT( complex *c, int opt ) {for( int i = 0;i < len;i ++ ) if( i < r[i] ) swap( c[i], c[r[i]] );for( int i = 1;i < len;i <<= 1 ) {complex omega( cos( pi / i ), sin( pi / i ) * opt );for( int j = 0;j < len;j += ( i << 1 ) ) {complex w( 1, 0 );for( int k = 0;k < i;k ++, w = w * omega ) {complex x = c[j + k], y = c[j + k + i] * w;c[j + k] = x + y, c[j + k + i] = x - y;}}}if( opt == -1 ) for( int i = 0;i < len;i ++ ) c[i].x /= len;
}#define inf 0x3f3f3f3f
int n, m, k;
int cnt[maxn];
char s[maxn], t[maxn];void work( char ch ) {for( int i = 0;i < len;i ++ ) f[i] = g[i] = complex( 0, 0 );int lst = -inf;for( int i = 0;i < n;i ++ ) {if( s[i] == ch ) lst = i;if( i - lst <= k ) f[i] = complex( 1, 0 );}lst = inf;for( int i = n - 1;~ i;i -- ) {if( s[i] == ch ) lst = i;if( lst - i <= k ) f[i] = complex( 1, 0 );}for( int i = 0;i < m;i ++ )if( t[i] == ch ) g[m - i - 1] = complex( 1, 0 );FFT( f, 1 ), FFT( g, 1 );for( int i = 0;i < len;i ++ ) f[i] = f[i] * g[i];FFT( f, -1 );for( int i = 0;i < len;i ++ ) cnt[i] += int( f[i].x + 0.5 );
}signed main() {scanf( "%d %d %d %s %s", &n, &m, &k, s, t );for( len = 1;len <= n + m;len <<= 1, l ++ );for( int i = 0;i < len;i ++ ) r[i] = r[i >> 1] >> 1 | ( i & 1 ) << l - 1;work( 'A' ), work( 'T' ), work( 'G' ), work( 'C' );int ans = 0;for( int i = 0;i < len;i ++ ) if( cnt[i] == m ) ans ++;printf( "%d\n", ans );return 0;
}
W:CodeForces 954I Yet Another String Matching Problem
洛谷链接
两个字符串比对,通过一些操作使得相应位是一样的。这种“比对”一般都是转化为图论,然后与连通块有关。
基础模型可以认为是两个序列,连边 Ai→BiA_i\rightarrow B_iAi→Bi 求连通块个数。
先只考虑两个串的距离如何计算?
连边 Si→TiS_i\rightarrow T_iSi→Ti (忽略原本就相同的位置)后,会形成若干的连通块。
一个联通块中只需要操作 ∣连通块∣−1|连通块|-1∣连通块∣−1 次就能将整个连通块变为相同的字符。
多个连通块之间是有向边的存在,是一个有向无环图,最后会归中于某个字符。
所以两个串的距离即为 (∑∣连通块∣)−1\Big(\sum{|连通块|}\Big)-1(∑∣连通块∣)−1。
具体实现可以用并查集维护字符集 {a,b,c,d,e,f}\{a,b,c,d,e,f\}{a,b,c,d,e,f}。每次合并原来不同的两个字符集的时候,答案 +1+1+1 即可。
现在考虑有多个 SSS 的子集 S′S'S′ 在与 TTT 适配。
设 ansi:ans_i:ansi: 考虑 S[i−∣T∣+1,i]S[i-|T|+1,i]S[i−∣T∣+1,i] 与 TTT 的适配距离答案。
同样地,分字符计算贡献再累和。
枚举 S′S'S′ 的字符 xxx,和 TTT 的字符 yyy 进行适配,x≠yx\neq yx=y。
我们需要判断是否有 Si′=x,Ti=yS'_i=x,T_i=ySi′=x,Ti=y 的位置,如果存在就意味着是需要操作的。
就执行并查集合并,答案 +1+1+1,如果之前的某些字符配对组合已经让 x,yx,yx,y 在同一个并查集,就不要加这个贡献。
每个 iii 不同结尾,适配距离也不同,对应的并查集也不同。所以每个位置都要自开一个独立的并查集。
本题的距离适配完全是适配的基础模型,翻转 TTT 再卷积为 fff,判断 fif_ifi 是否为 111,在 fif_ifi 的并查集里面进行合并即可。
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
#define maxn 1000000
int len, l, n, m;
int r[maxn], ans[maxn];
char s[maxn], t[maxn];struct complex {double x, i;complex(){}complex( double X, double I ) { x = X; i = I; }
}f[maxn], g[maxn];
double pi = acos( -1.0 );
complex operator + ( complex u, complex v ) {return complex( u.x + v.x, u.i + v.i );
}
complex operator - ( complex u, complex v ) {return complex( u.x - v.x, u.i - v.i );
}
complex operator * ( complex u, complex v ) {return complex( u.x * v.x - u.i * v.i, u.x * v.i + v.x * u.i );
}void FFT( complex *c, int opt ) {for( int i = 0;i < len;i ++ ) if( i < r[i] ) swap( c[i], c[r[i]] );for( int i = 1;i < len;i <<= 1 ) {complex omega( cos( pi / i ), opt * sin( pi / i ) );for( int j = 0;j < len;j += ( i << 1 ) ) {complex w( 1, 0 );for( int k = 0;k < i;k ++, w = w * omega ) {complex x = c[j + k], y = c[j + k + i] * w;c[j + k] = x + y, c[j + k + i] = x - y;}}}if( opt == -1 ) for( int i = 0;i < len;i ++ ) c[i].x /= len;
}struct DSU {int f[10];void init() { for( int i = 0;i < 6;i ++ ) f[i] = i; }int find( int x ) { return f[x] == x ? x : f[x] = find( f[x] ); }void merge( int x, int y ) { x = find( x ), y = find( y ), f[y] = x; }
}dsu[maxn];void work( int x, int y ) {for( int i = 0;i < len;i ++ ) f[i] = g[i] = complex( 0, 0 );for( int i = 0;i < n;i ++ ) f[i].x = ( s[i] - 'a' == x );for( int i = 0;i < m;i ++ ) g[m - i - 1].x = ( t[i] - 'a' == y );FFT( f, 1 ), FFT( g, 1 );for( int i = 0;i < len;i ++ ) f[i] = f[i] * g[i];FFT( f, -1 );for( int i = 0;i < n;i ++ ) {int c = int( f[i].x + 0.5 );if( ! c ) continue;if( dsu[i].find( x ) ^ dsu[i].find( y ) ) ans[i] ++, dsu[i].merge( x, y );}
}int main() {scanf( "%s %s", s, t );n = strlen( s ), m = strlen( t );for( len = 1;len <= n + m;len <<= 1, l ++ );for( int i = 0;i < len;i ++ ) r[i] = r[i >> 1] >> 1 | ( i & 1 ) << l - 1;for( int i = 0;i < n;i ++ ) for( int j = 0;j < 6;j ++ ) dsu[i].init();for( int i = 0;i < 6;i ++ ) for( int j = 0;j < 6;j ++ ) if( i ^ j ) work( i, j );for( int i = m - 1;i < n;i ++ ) printf( "%d\n", ans[i] );return 0;
}