当前位置:网站首页>[cf1054h] epic Revolution -- number theory, convolution, arbitrary modulus NTT
[cf1054h] epic Revolution -- number theory, convolution, arbitrary modulus NTT
2022-07-29 07:00:00 【Ouye xjx】
CF1054H Epic Convolution
Answer key
This question is called Epic Convolution, Which translates as “ Epic convolution ”( Also called “ Shit shit convolution ”)
Because the modulus is very small and a prime number , So by Euler's theorem , We can turn the problem into a quest
f x = ∑ i = 0 n − 1 ∑ j = 0 m − 1 A i B j [ i 2 j 3 m o d 490018 = x ] f_x=\sum_{i=0}^{n-1}\sum_{j=0}^{m-1}A_iB_j[i^2j^3\bmod 490018=x] fx=i=0∑n−1j=0∑m−1AiBj[i2j3mod490018=x] Then in order to solve more efficiently , We need to consider transforming the above multiplicative convolution into additive convolution .
can 490018 490018 490018 Obviously, there is no original root ? How to convert ?
At this time, the idea of dividing the range into Multiplication rings is obviously unrealistic ( That's why I'm a mathematical shit ), We need more abstract ideas . Consider putting 490018 490018 490018 Make a unique decomposition :
490018 = 2 × 491 × 499 490018=2\times 491 \times 499 490018=2×491×499 And when we get i 2 j 3 i^2j^3 i2j3 modulus 2 , 491 , 499 2,491,499 2,491,499 After the value of , It can be used CRT combined , So we might as well change the formula , Change to seek
f x , y , z = ∑ i = 0 n − 1 ∑ j = 0 m − 1 A i B j [ i 2 j 3 m o d 2 = x ] [ i 2 j 3 m o d 491 = y ] [ i 2 j 3 m o d 499 = z ] f_{x,y,z}=\sum_{i=0}^{n-1}\sum_{j=0}^{m-1}A_iB_j[i^2j^3\bmod 2=x][i^2j^3\bmod 491=y][i^2j^3\bmod 499=z] fx,y,z=i=0∑n−1j=0∑m−1AiBj[i2j3mod2=x][i2j3mod491=y][i2j3mod499=z]
The rest is simple , We just need to convolute each dimension in turn , The first dimension can be violently convoluted , The latter two dimensions can be transformed into additive convolution by taking logarithm of the original root .
Notice that 0 Cases where logarithms cannot be taken , stay 0 Place the convolution of violence .
Code
What I write is the convolution of violence classification discussion , In fact, a little bit of exclusion can significantly reduce the number of operations .
#include<bits/stdc++.h>//JZM yyds!!
#define ll long long
#define lll __int128
#define uns unsigned
#define fi first
#define se second
#define IF (it->fi)
#define IS (it->se)
#define lowbit(x) ((x)&-(x))
#define END putchar('\n')
#define inline jzmyyds
using namespace std;
const int MAXN=5e5+5;
const ll INF=1e17;
ll read(){
ll x=0;bool f=1;char s=getchar();
while((s<'0'||s>'9')&&s>0){
if(s=='-')f^=1;s=getchar();}
while(s>='0'&&s<='9')x=(x<<1)+(x<<3)+(s^48),s=getchar();
return f?x:-x;
}
int ptf[50],lpt;
void print(ll x,char c='\n'){
if(x<0)putchar('-'),x=-x;
ptf[lpt=1]=x%10;
while(x>9)x/=10,ptf[++lpt]=x%10;
while(lpt)putchar(ptf[lpt--]^48);
if(c>0)putchar(c);
}
//490018=2*491*499,g_491=2,g_499=7
const ll MOD=490019,NM=1000000000245104641ll;// One pass large modulus NTT
const lll E=1;
ll ksml(ll a,ll b,ll mo=NM){
ll res=1;
for(;b;b>>=1,a=E*a*a%mo)if(b&1)res=E*res*a%mo;
return res;
}
ll ksm(ll a,ll b,ll mo=MOD){
ll res=1;
for(;b;b>>=1,a=a*a%mo)if(b&1)res=res*a%mo;
return res;
}
#define g 7ll
int rev[1025];
int initrev(){
for(int i=0;i<(1<<10);i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<9);
return 114514;
}
ll omg[1025],cbddl=initrev(),IVN=ksml(1024,NM-2);
void NTT(ll*a,int inv,int n=1024){
for(int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
ll x,y,tmp=233,mi=(NM-1)>>1;omg[0]=1;
for(int m=1;m<n;m<<=1,mi>>=1){
if(m>1)tmp=ksml(g,inv>0?mi:NM-1-mi);
for(int i=1;i<m;i++)omg[i]=E*omg[i-1]*tmp%NM;
for(int i=0,om=0;i<n;i+=(m<<1),om=0)
for(int j=i;j<i+m;j++,om++){
x=a[j],y=E*a[j+m]*omg[om]%NM,a[j]=x+y,a[j+m]=x-y;
if(a[j]>=NM)a[j]-=NM;
if(a[j+m]<0)a[j+m]+=NM;
}
}if(inv<0)for(int i=0;i<n;i++)a[i]=E*a[i]*IVN%NM;
}
#undef g
ll f[1025][1025],g[1025][1025];
void polymul(){
for(int i=0;i<500;i++)NTT(f[i],1),NTT(g[i],1);
for(int i=1;i<1024;i++)
for(int j=0;j<i;j++)swap(f[i][j],f[j][i]),swap(g[i][j],g[j][i]);
for(int i=0;i<1024;i++){
NTT(f[i],1),NTT(g[i],1);
for(int j=0;j<1024;j++)f[i][j]=E*f[i][j]*g[i][j]%NM;
NTT(f[i],-1);
}
for(int i=1;i<1024;i++)for(int j=0;j<i;j++)swap(f[i][j],f[j][i]);
for(int i=0;i<1000;i++)NTT(f[i],-1);
}
int ex1[810],lg1[810],ex2[810],lg2[810];
int initexp(){
ex1[0]=ex2[0]=1,lg1[1]=lg2[1]=0;
for(int i=1;i<490;i++)ex1[i]=(ex1[i-1]<<1)%491,lg1[ex1[i]]=i;
for(int i=1;i<498;i++)ex2[i]=ex2[i-1]*7%499,lg2[ex2[i]]=i;
return 1919810;
}
int sh_t=initexp(),n,m;
ll a1[499],b1[499],ans,mi[MAXN],A[MAXN];
void solve(ll(*a)[499],ll(*b)[499],ll(*c)[499],bool CL){
if(CL){
memset(f,0,sizeof(f));
memset(g,0,sizeof(g));
memset(a1,0,sizeof(a1));
memset(b1,0,sizeof(b1));
}
for(int i=1;i<491;i++)
for(int j=1;j<499;j++)
(f[lg1[i]][lg2[j]]+=a[i][j])%=MOD,(g[lg1[i]][lg2[j]]+=b[i][j])%=MOD;
for(int i=1;i<491;i++)
for(int j=1;j<499;j++)(a1[i]+=a[i][j])%=MOD;
for(int i=1;i<491;i++)
for(int j=0;j<499;j++)(b1[i]+=b[i][j])%=MOD;
polymul();
for(int i=0,mi=0;i<1000;i++,mi=i%490)
for(int j=0;j<1000;j++)
(c[ex1[mi]][ex2[j%498]]+=f[i][j])%=MOD;
for(int i=1;i<491;i++)
for(int j=1,to=i;j<491;j++,to=(to+i)%491)
(c[to][0]+=a[i][0]*b1[j]+a1[i]*b[j][0])%=MOD;
memset(a1,0,sizeof(a1));
memset(b1,0,sizeof(b1));
for(int i=1;i<491;i++)
for(int j=0;j<499;j++)(a1[j]+=a[i][j])%=MOD;
for(int i=0;i<491;i++)
for(int j=0;j<499;j++)(b1[j]+=b[i][j])%=MOD;
for(int i=0;i<499;i++)
for(int j=0,to=0;j<499;j++,to=(to+i)%499)
(c[0][to]+=a[0][i]*b1[j]+a1[i]*b[0][j])%=MOD;
}
int crt(int a,int b,int c){
static const ll iv1=492>>1,iv2=ksm(491<<1,499-2,499);
int k1=(b-a+491)*iv1%491,k2=(c-a-k1-k1+4990)*iv2%499;
return (a+(k1<<1)+(k2*491<<1))%(MOD-1);
}
ll a[2][491][499],b[2][491][499],b_[491][499],as[2][491][499],c;
const bool KFC=1;
int main()
{
n=read(),m=read(),c=read(),mi[0]=1;
for(int i=1;i<MOD;i++)mi[i]=mi[i-1]*c%MOD;
if(1ll*n*m<=1000000ll&&KFC){
for(int i=0;i<n;i++)A[i]=read();
for(int j=0;j<m;j++){
ll B=read();
for(int i=0;i<n;i++)
(ans+=A[i]*B*mi[1ll*i*i*j%(MOD-1)*j*j%(MOD-1)])%=MOD;
}return print(ans),0;
}
for(int i=0,k=0;i<n;i++,k=1ll*i*i%(MOD-1))
(a[k&1][k%491][k%499]+=read())%=MOD;
for(int i=0,k=0,d;i<m;i++,k=1ll*i*i*i%(MOD-1))
d=read(),(b[k&1][k%491][k%499]+=d)%=MOD,(b_[k%491][k%499]+=d)%=MOD;
solve(a[1],b[1],as[1],0);
solve(a[0],b_,as[0],1);
solve(a[1],b[0],as[0],1);
for(int i=0;i<2;i++)for(int j=0;j<491;j++)
for(int k=0;k<499;k++)(ans+=as[i][j][k]*mi[crt(i,j,k)])%=MOD;
print(ans);
return 0;
}
边栏推荐
- Unity exploration plot access design analysis & process + code specific implementation
- 5g service interface and reference point
- SSH免密登录-两台虚拟机建立免密通道 双向信任
- 网上传说软件测试培训真的那么黑心吗?都是骗局?
- Etcd principle
- Unity探索地块通路设计分析 & 流程+代码具体实现
- 模拟卷Leetcode【普通】222. 完全二叉树的节点个数
- Talk about tcp/ip protocol? And the role of each layer?
- 实战!聊聊如何解决MySQL深分页问题
- 数据库持久化+JDBC数据库连接
猜你喜欢

MySQL:当你CRUD时BufferPool中发生了什么?十张图就能说清楚

LDAP brief description and unified authentication description

吴恩达老师机器学习课程笔记 02 单变量线性回归

10 frequently asked JVM questions in interviews

SDN topology discovery principle

Why does 5g N2 interface control plane use SCTP protocol?

IO stream - file - properties

Sword finger offer II 115: reconstruction sequence

分享一些你代码更好的小建议,流畅编码提搞效率

微信小程序的反编译
随机推荐
【干货备忘】50种Matplotlib科研论文绘图合集,含代码实现
Pytorch多GPU条件下DDP集群分布训练实现(简述-从无到有)
Leetcode-592: fraction addition and subtraction
Teacher Wu Enda's machine learning course notes 02 univariate linear regression
mysql可以定时导出表格吗?
关于SQL Server语句入门级应用阶段性学习——找工作必备(一)
数据库多表查询 联合查询 增删改查
The core of openresty and cosocket
【论文阅读 | cryoET】Gum-Net:快速准确的3D Subtomo图像对齐和平均的无监督几何匹配
Etcd principle
竣达技术 | 适用于”日月元”品牌UPS微信云监控卡
【笔记】The art of research(明白问题的重要性)
没那么简单的单例模式
游戏资产的革命
Simulation volume leetcode [ordinary] 172. Zero after factorial
数据库持久化+JDBC数据库连接
Security in quantum machine learning
DBAsql面试题
Teacher Wu Enda machine learning course notes 01 introduction
Cesium reflection