码迷,mamicode.com
首页 > 其他好文 > 详细

COGS2294 释迦

时间:2017-02-25 16:56:13      阅读:257      评论:0      收藏:0      [点我收藏+]

标签:namespace   out   级别   names   alt   href   display   注意   const   

传送门

就是传说中的任意模数卷积嘛……有三模数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(\mod m_1)\\Ans\equiv a_2(\mod m_2)\\Ans\equiv a_3(\mod m_3)\end{equation}

先用中国剩余定理合并前两个同余式,得到

\begin{equation}Ans\equiv A(\mod M)\\Ans\equiv a_3(\mod m_3)\end{equation}

不妨设

\begin{equation}Ans=kM+A=xm_3+a_3\end{equation}

我们可以在$\mod m_3$的意义下求解k的值,那么有

\begin{equation}kM\equiv a_3-A(\mod m_3)\end{equation}

(因为是在$\mod m_3$的意义下,所以$xm_3$被消掉了)

也就是说

\begin{equation}k\equiv (a_3-A)M^{-1}(\mod m_3)\end{equation}

求出$k$之后代入$Ans=kM+A$,这次只要在$\mod 23333333$的意义下算出$Ans$的值即可。

技术分享
 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 using namespace std;
 5 const int maxn=262200,m1=998244353,m2=1004535809,m3=469762049,g=3,Mod=23333333;
 6 const long long M=(long long)m1*m2;
 7 void NTT(int*,int,int,int);
 8 int China(int,int,int);
 9 int qpow(int,int,int);
10 long long mul(long long,long long,long long);
11 int n,N=1,A[maxn],B[maxn],C[maxn],D[maxn],a[3][maxn];
12 int main(){
13     freopen("annona_squamosa.in","r",stdin);
14     freopen("annona_squamosa.out","w",stdout);
15     scanf("%d",&n);
16     while(N<(n<<1))N<<=1;
17     for(int i=0;i<n;i++)scanf("%d",&A[i]);
18     for(int i=0;i<n;i++)scanf("%d",&B[i]);
19     copy(A,A+N,C);
20     copy(B,B+N,D);
21     NTT(C,N,1,m1);
22     NTT(D,N,1,m1);
23     for(int i=0;i<N;i++)a[0][i]=(long long)C[i]*D[i]%m1;
24     NTT(a[0],N,-1,m1);
25     copy(A,A+N,C);
26     copy(B,B+N,D);
27     NTT(C,N,1,m2);
28     NTT(D,N,1,m2);
29     for(int i=0;i<N;i++)a[1][i]=(long long)C[i]*D[i]%m2;
30     NTT(a[1],N,-1,m2);
31     copy(A,A+N,C);
32     copy(B,B+N,D);
33     NTT(C,N,1,m3);
34     NTT(D,N,1,m3);
35     for(int i=0;i<N;i++)a[2][i]=(long long)C[i]*D[i]%m3;
36     NTT(a[2],N,-1,m3);
37     for(int i=0;i<n;i++)printf("%d\n",China(a[0][i],a[1][i],a[2][i]));
38     return 0;
39 }
40 void NTT(int *A,int n,int tp,int p){
41     for(int i=0;i<n;i++)A[i]%=p;
42     for(int i=1,j=0,k;i<n-1;i++){
43         k=n;
44         do j^=(k>>=1);while(j<k);
45         if(i<j)swap(A[i],A[j]);
46     }
47     for(int k=2;k<=n;k<<=1){
48         int wn=qpow(g,(tp>0?(p-1)/k:(p-1)/k*(long long)(p-2)%(p-1)),p);
49         for(int i=0;i<n;i+=k){
50             int w=1;
51             for(int j=0;j<(k>>1);j++,w=(long long)w*wn%p){
52                 int a=A[i+j],b=(long long)w*A[i+j+(k>>1)]%p;
53                 A[i+j]=(a+b)%p;
54                 A[i+j+(k>>1)]=(a-b+p)%p;
55             }
56         }
57     }
58     if(tp<0){
59         int inv=qpow(n,p-2,p);
60         for(int i=0;i<n;i++)A[i]=(long long)A[i]*inv%p;
61     }
62 }
63 int China(int a1,int a2,int a3){
64     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;
65     return ((k%Mod)*(M%Mod)%Mod+A%Mod)%Mod;
66 }
67 int qpow(int a,int b,int p){
68     int ans=1;
69     for(;b;b>>=1,a=(long long)a*a%p)if(b&1)ans=(long long)ans*a%p;
70     return ans;
71 }
72 long long mul(long long a,long long b,long long p){
73     a%=p;b%=p;
74     return ((a*b-(long long)((long long)((long double)a/p*b+1e-3)*p))%p+p)%p;
75 }
View Code

 

COGS2294 释迦

标签:namespace   out   级别   names   alt   href   display   注意   const   

原文地址:http://www.cnblogs.com/hzoier/p/6441967.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!