COGS2294 释迦
就是传说中的任意模数卷积嘛……有三模数NTT和拆系数FFT等做法,我比较懒不想动脑子,就用了三模数NTT的做法……
卷积之后每个数可以达到$10^{23}$左右的级别,直接long double或者__float128都会炸精度(而且__float128炸得更惨……好像是转换的时候掉精度太多……)。而这个模数又不能NTT(首先这就不是个质数……),因此直接搞是行不通的。
我们可以选三个满足NTT性质并且乘起来$>10^{23}$的模数分别NTT,最后中国剩余定理合并。但注意到$10^{23}>2^{64}$,因此直接合并会炸long long,所以我们就需要一些tricky的办法来合并。
我们得到的是这样的三个同余式:
\begin{equation}Ans\equiv a_1\pmod{m_1}\\Ans\equiv a_2\pmod{m_2}\\Ans\equiv a_3\pmod{m_3}\end{equation}
先用中国剩余定理合并前两个同余式,得到
\begin{equation}Ans\equiv A{\pmod M}\\Ans\equiv a_3\pmod{m_3}\end{equation}
不妨设
\begin{equation}Ans=kM+A=xm_3+a_3\end{equation}
我们可以在$\bmod m_3$的意义下求解k的值,那么有
\begin{equation}kM\equiv a_3-A\pmod{m_3}\end{equation}
(因为是在$\bmod m_3$的意义下,所以$xm_3$被消掉了)
也就是说
\begin{equation}k\equiv (a_3-A)M^{-1}\pmod{m_3}\end{equation}
求出$k$之后代入$Ans=kM+A$,这次只要在$\bmod 23333333$的意义下算出$Ans$的值即可。
#include<cstdio> #include<cstring> #include<algorithm> using namespace std; const int maxn=262200,m1=998244353,m2=1004535809,m3=469762049,g=3,Mod=23333333; const long long M=(long long)m1*m2; void NTT(int*,int,int,int); int China(int,int,int); int qpow(int,int,int); long long mul(long long,long long,long long); int n,N=1,A[maxn],B[maxn],C[maxn],D[maxn],a[3][maxn]; int main(){ freopen("annona_squamosa.in","r",stdin); freopen("annona_squamosa.out","w",stdout); scanf("%d",&n); while(N<(n<<1))N<<=1; for(int i=0;i<n;i++)scanf("%d",&A[i]); for(int i=0;i<n;i++)scanf("%d",&B[i]); copy(A,A+N,C); copy(B,B+N,D); NTT(C,N,1,m1); NTT(D,N,1,m1); for(int i=0;i<N;i++)a[0][i]=(long long)C[i]*D[i]%m1; NTT(a[0],N,-1,m1); copy(A,A+N,C); copy(B,B+N,D); NTT(C,N,1,m2); NTT(D,N,1,m2); for(int i=0;i<N;i++)a[1][i]=(long long)C[i]*D[i]%m2; NTT(a[1],N,-1,m2); copy(A,A+N,C); copy(B,B+N,D); NTT(C,N,1,m3); NTT(D,N,1,m3); for(int i=0;i<N;i++)a[2][i]=(long long)C[i]*D[i]%m3; NTT(a[2],N,-1,m3); for(int i=0;i<n;i++)printf("%d\n",China(a[0][i],a[1][i],a[2][i])); return 0; } void NTT(int *A,int n,int tp,int p){ for(int i=0;i<n;i++)A[i]%=p; for(int i=1,j=0,k;i<n-1;i++){ k=n; do j^=(k>>=1);while(j<k); if(i<j)swap(A[i],A[j]); } for(int k=2;k<=n;k<<=1){ int wn=qpow(g,(tp>0?(p-1)/k:(p-1)/k*(long long)(p-2)%(p-1)),p); for(int i=0;i<n;i+=k){ int w=1; for(int j=0;j<(k>>1);j++,w=(long long)w*wn%p){ int a=A[i+j],b=(long long)w*A[i+j+(k>>1)]%p; A[i+j]=(a+b)%p; A[i+j+(k>>1)]=(a-b+p)%p; } } } if(tp<0){ int inv=qpow(n,p-2,p); for(int i=0;i<n;i++)A[i]=(long long)A[i]*inv%p; } } int China(int a1,int a2,int a3){ long long A=(mul((long long)a1*m2%M,qpow(m2%m1,m1-2,m1),M)+mul((long long)a2*m1%M,qpow(m1%m2,m2-2,m2),M))%M,k=((a3-A)%m3+m3)%m3*qpow(M%m3,m3-2,m3)%m3; return ((k%Mod)*(M%Mod)%Mod+A%Mod)%Mod; } int qpow(int a,int b,int p){ int ans=1; for(;b;b>>=1,a=(long long)a*a%p)if(b&1)ans=(long long)ans*a%p; return ans; } long long mul(long long a,long long b,long long p){ a%=p;b%=p; return ((a*b-(long long)((long long)((long double)a/p*b+1e-3)*p))%p+p)%p; }View Code
- [任意模数NTT 三模数NTT] COGS 2294 [HZOI 2015] 释迦
- COGS 2294. [HZOI 2015] 释迦 (FFT mod any prime)
- COGS 614 游历校园
- cogs896. 圈奶牛
- 【COGS】256 [POI2001] 金矿 线段树
- COGS 1190 最大数
- cogs爱争吵的猴子 题解
- [cogs] 1767 [NOI2014]随机数生成器
- 【线段树】【NOI 1999】【cogs 284】内存分配
- NOIP 2012 国王游戏 贪心 高精度 (COGS 1263)
- 数列 COGS1048:[Citric S2] 一道防AK好题
- cogs 547:[HAOI2011] 防线修建
- 图论(对偶图):COGS 470. [NOI2010]海拔
- 数据结构(主席树):COGS 2211. 谈笑风生
- cogs 自己出的题目 题解报告
- COGS 137. [USACO Feb08] 连线游戏
- 字符串(后缀自动机):COGS 2399. 循环同构
- 【BIT套主席树】COGS257-动态排名系统
- COGS 2060
- cogs 495 窗口 单调队列第一题