Loj #6019. 「from CommonAnts」寻找 LCM

给个链接:https://loj.ac/problem/6019

 

还是一道扩展卢卡斯+中国剩余定理。。。。就当练练手

但是这题怎么这么卡常呢????!!!!!

在LOJ上死也过不去 (为什么要加那么多不自然的优化才能过2333考试的时候谁能想到)

 

然后在洛谷上测了一下,果然被卡了四倍多常数哎。。。。。

 

 

#include<bits/stdc++.h>
#define maxn 1000005
#define ll long long
using namespace std;
int d[30],c[30],D[30];
int ans[30],n,m,P,mo,MOD;
int a[maxn],b[maxn],num;
int jc[maxn],v[maxn],dfn;

inline int add(int x,int y,const int ha){
	x+=y;
	if(x>=ha) return x-ha;
	else return x;
}

inline int ksm(int x,int y,const int ha){
	int an=1;
	for(;y;y>>=1,x=x*(ll)x%ha) if(y&1) an=an*(ll)x%ha;
	return an;
}

inline int getni(int x){
	return ksm(x,D[mo]/d[mo]*(d[mo]-1)-1,D[mo]);
}

struct node{
	int val,tmp;
	node operator *(const node &U)const{
		return (node){val*(ll)U.val%D[mo],tmp+U.tmp};
	}
	node operator /(const node &U)const{
		return (node){val*(ll)getni(U.val)%D[mo],tmp-U.tmp};
	}
}ANS[maxn];

inline void dvd(){
	MOD=P;
	for(int i=2;i*(ll)i<=P;i++) if(!(P%i)){
		d[++num]=i,D[num]=1;
		while(!(P%i)) P/=i,D[num]*=i,c[num]++;
		if(P==1) break;
	}
	if(P!=1) d[++num]=D[num]=P,c[num]=1;
}

inline node getjc(int x){
	if(v[x]==dfn) return ANS[x];
	
	node now={1,0};
	if(x>=d[mo]) now=now*getjc(x/d[mo]),now.tmp+=x/d[mo];
	if(x>=D[mo]) now=now*(node){ksm(jc[D[mo]-1],x/D[mo],D[mo]),0};
	now=now*(node){jc[x%D[mo]],0};
	
	ANS[x]=now,v[x]=dfn;
	return now;
}

inline node getC(int x,int y){
	return getjc(x)/getjc(y)/getjc(x-y);
}

inline void solve(int x){
	mo=x,jc[0]=1,dfn++;
	int tp=min(1000000,D[x]-1);
	for(int i=1;i<=tp;i++){
		jc[i]=jc[i-1];
		if(i%d[x]) jc[i]=jc[i]*(ll)i%D[x];
	}
	
	ans[x]=1;
	node now;
	for(int i=1;i<=n;i++){
		now=getC(b[i],a[i]);
		ans[x]=ans[x]*(ll)now.val%D[x]*(ll)ksm(d[x],now.tmp,D[x])%D[x];
	}
}

inline int CRT(){
	int an=0;
	for(int i=1;i<=num;i++){
		mo=i;
		an=add(an,getni(MOD/D[i])*(ll)(MOD/D[i])%MOD*(ll)ans[i]%MOD,MOD);
	}
	return an;
}

int main(){
	scanf("%d%d",&n,&P);
	for(int i=1;i<=n;i++) scanf("%d",a+i);
	for(int i=1;i<=n;i++) scanf("%d",b+i);
	dvd();
	for(int i=1;i<=num;i++) solve(i);
	printf("%d\n",CRT());
	return 0;
}

  

posted @ 2018-02-22 13:00  蒟蒻JHY  阅读(1122)  评论(0编辑  收藏  举报