题解 LOJ3050 「十二省联考 2019」骗分过样例

CASE (1sim 3)

(n)组测试数据,每次输入一个数(x),求(19^x)

测试点(1)(x=0,1,dots n-1),可以直接递推。

测试点(2)要开long long并用快速幂。

测试点(3)(x)超出了long long范围。根据欧拉定理,当(a,p)互质时,(a^bequiv a^{bmodvarphi(p)}pmod p)。因为模数(p=998244353)是一个质数,所以(varphi(p)=p-1=998244352)。将数字一位一位读进来然后对(998244352)取模即可。

代码片段:

typedef long long ll;
ll pow_mod(...){
	//快速幂略
}
ll read_mod(){
	ll x=0;
	char ch=getchar();
	while(!isdigit(ch))ch=getchar();
	while(isdigit(ch))x=(x*10+ch-'0')%Math::PHI_MOD,ch=getchar();
	return x;
}
void solve1(ll x=998244353){
	Init(x);//MOD=x
	int n=read();while(n--){
		ll x=read_mod();
		printf("%lld
",pow_mod(19,x));
	}
}

CASE (4sim 5)

测试点(4),功能编号是1?,容易猜到意思是说模数不确定。于是考虑枚举模数。先找出输出文件中最大的数(mx),则模数一定大于(mx)。我们从(mx+1)开始枚举模数。容易找到模数为(1145141)

找模数代码:

注意这里我们只判断了模数是质数的情况,因为我们相信出题人足够良心

ll read_mod(string s){
	ll x=0;
	for(int i=0;i<(int)s.length();++i)x=(x*10+s[i]-'0')%Math::PHI_MOD;
	return x;
}
bool test_mod(ll m,vector<pair<string,ll> >&vec){
	if(!is_prime(m))return false;
	Init(m);
	for(int i=0;i<(int)vec.size();++i){
		ll x=read_mod(vec[i].fst);
		if(pow_mod(19,x)!=vec[i].scd)return false;
	}
	return true;
}
void get_mod(){
	freopen("software4.in","r",stdin);
	int n;
	string str;cin>>str>>n;
	for(int i=1;i<=3;++i)cin>>str;
	vector<pair<string,ll> >vec;
	for(int i=0;i<5;++i){
		cin>>str;
		vec.pb(mk(str,233));
	}
	freopen("software4.ans","r",stdin);
	for(int i=1;i<=3;++i)cin>>str;
	ll mx=0;
	for(int i=0;i<n;++i){
		ll x;cin>>x;
		if(i<(int)vec.size())vec[i].scd=x;
		mx=max(mx,x);
	}
//	for(int i=0;i<(int)vec.size();++i){
//		cerr<<vec[i].fst<<" "<<vec[i].scd<<endl;
//	}
	cerr<<"max = "<<mx<<endl;
	for(ll i=mx+1;i<=LLONG_MAX;++i){
		if(test_mod(i,vec)){
			cerr<<i<<endl;
			return;
		}
	}
}

找到模数后:

void solve2(){
	solve1(1145141);
}

测试点(5),发现答案都很大,说明模数也很大,用刚刚的方法大概是别想在考试结束前找到模数了。

考虑更巧妙的数论技巧。我们在输入文件中找两个接近的数(x,y),满足(x<y)(ans_x>ans_y)。假设模数为(p),则可以知道(ans_xcdot19^{y-x}equiv ans_ypmod p)

所以模数(p)就一定是(ans_xcdot 19^{y-x}-ans_y)的约数。注意这个数本身很大,要用__int128存。暴力枚举其所有约数复杂度是(O(sqrt{n}))的,依然无法接受。于是做如下剪枝:

  • 从输出文件中找出最大的数(mx),模数一定大于(mx)

  • 因为我们相信出题人足够良心,所以模数应该是质数。

  • 因为我们相信出题人足够良心,所以模数应该在long long范围内。

判断大质数的方法会在CASE 8~10部分详细介绍。

用如下的程序可以很快地求出,模数是(5211600617818708273)

pair<int,ll>data[10004];
void print128(__int128 x){
	const ll base=1e18;
	if(x<base)cout<<(ull)x<<endl;
	else{
		cout<<(ll)(x/base);
		ll t=x%base;int len=0;
		while(t)len++,t/=10;
		for(int i=1;i<=18-len;++i)cout<<"0";
		if(x%base)cout<<(ll)(x%base)<<endl;
	}
}
void get_big_mod(){
	freopen("software5.in","r",stdin);
	string s;cin>>s;
	int n;cin>>n;
	for(int i=1;i<=n;++i)cin>>data[i].fst;
	freopen("software5.ans","r",stdin);
	ll mx=0;
	for(int i=1;i<=n;++i){cin>>data[i].scd;mx=max(mx,data[i].scd);}
	sort(data+1,data+n+1);
	int x=0,y=0;
	for(int i=1;i+1<=n;++i){
		for(int j=i+1;j<=n;++j){
			if(data[j].fst==data[i].fst)continue;
			if(data[j].scd>data[i].scd)continue;
			// i.f<j.f && i.s>j.s
			if(!x||data[j].fst-data[i].fst<data[y].fst-data[x].fst)x=i,y=j;
		}
	}
	cerr<<data[x].fst<<" "<<data[x].scd<<endl;
	cerr<<data[y].fst<<" "<<data[y].scd<<endl;
	cerr<<"--------------------"<<endl;
	__int128 X=(__int128)data[x].scd*(__int128)19*(__int128)19-(__int128)data[y].scd;
	print128(X);
	cerr<<"--------------------"<<endl;
	vector<pair<string,ll> >vec;
	for(int i=n;i>=n-5+1;--i){
		string str="";
		ull x=data[i].fst;
		while(x)str=(char)('0'+x%10)+str,x/=10;
		vec.pb(mk(str,data[i].scd));
	}
	for(__int128 i=2;i*i<=X;++i){
		if(X%i!=0)continue;
		if(i>mx&&i<LLONG_MAX){
			if(test_mod(i,vec)){
				cerr<<(ll)i<<endl;
			}
		}
		if(X/i>mx&&X/i<LLONG_MAX){
			if(test_mod(X/i,vec)){
				cerr<<(ll)(X/i)<<endl;
			}
		}
	}
}

找到模数后:

void solve3(){
	solve1(5211600617818708273LL);
}

CASE (6sim 7)

功能编号是wa,注意到前几个数和第一个测试点相同,后面开始出现负数。结合题目下方关于自然溢出的提示,不难猜到这个两个测试点应该是(19^x)发生了自然溢出

注意,自然溢出的过程非常奇特,不能用快速幂,否则就溢不出我们想要的效果。最保险的方法是按照题目的提示做加法。

对于测试点(7),输入的数值都非常大,直接做加法递推显然是不合适的。打表可以发现,(19^x)的自然溢出有循环节。这个循环节开始于(55245),长度是(45699)

找循环节的代码:

//#include <windows.h>
map<int,int>mp;
void get_repetend(){
	//找自然溢出循环节
	Init(998244353);
	int x=1;
	for(int i=0;i<10000000;++i){
		if(mp.count(x)==0)mp[x]=i;
		else{
			//cerr<<mp[x]<<" "<<i<<endl;//Sleep(100);
			cerr<<"start: "<<mp[x]<<endl<<"term: "<<i-mp[x]<<endl;
			return;
		}
		int new_x=x;
		for(int j=1;j<=18;++j)new_x=((int)((uint)new_x+(uint)x));
		x=new_x%Math::MOD;
	}
}

找出循环节之后的实现:

const int ST=55245,LEN=45699;
int ans[ST+LEN];
void solve4(){
	Init(998244353);//MOD=998244353
	int x=1;
	for(int i=0;i<ST+LEN;++i){
		ans[i]=x;
		int new_x=x;
		for(int j=1;j<=18;++j)new_x=((int)((uint)new_x+(uint)x));
		x=new_x%Math::MOD;
	}
	int n=read();
	for(int i=1;i<=n;++i){
		ll x=readll();
		if(x<ST+LEN)printf("%d
",ans[x]);
		else printf("%d
",ans[(x-ST)%LEN+ST]);
	}
}

CASE (8sim 10)

首先观察到,输入的每一行,对应地输出了一个字符串。设输入的数为(x,y),则输出字符串的长度为(y-x+1)。不难猜到输出的是([x,y])这一区间内每个数的信息。由功能编号p可以联想到prime number,质数。由此发现,输出字符串的每一位,表示([x,y])区间里的每个数是否是质数,若为质数则该位为p,否则为.

测试点(8)可以使用线性筛法。

对于测试点(9,10),这里介绍一种快速的素数判定算法:米勒-罗宾法。

前文我们已经介绍了欧拉定理,对于它的一种特殊情况:模数为质数时,我们得到费马小定理

(p)为质数,(forall a<p),有:

[a^{p-1}equiv 1pmod p ]

我们选取一些(a),然后把(p)带入,进行判定。若对于其中某个(a_i)(p)不满足费马小定理,则可以确定(p)不是质数。

但这个条件是必要不充分条件。也就是说,存在一类伪素数,它们虽然不是素数,但满足费马小定理。

为避免这种情况,我们的算法要进行二次探测。它基于如下定理:

对一个素数(p),若有(x^2equiv 1pmod p),则(xequiv 1pmod p)(xequiv p-1pmod p)

证明:移项得(x^2-1equiv 0pmod p)。用平方差公式展开得((x+1)(x-1)equiv 0pmod p)。所以(xequiv 1pmod p)(xequiv p-1pmod p)

很遗憾,这个定理同样是必要不充分条件。但我们可以利用它来加强我们的条件,提高素数判定的准确性。

请注意,接下来我们讨论的是一个递归的过程,在递归过程中,指数(p)是在不断变化的,而模数(P)是始终不变的。初始时,(P=p),都是我们要判定的那个数。

首先,经过费马小定理的验证,我们已经确定(a^{p-1}equiv 1pmod P)。如果(p-1)是偶数,则有(left(a^{frac{p-1}{2}} ight)^2equiv 1pmod P)。那么可以将(a^{frac{p-1}{2}})视作“二次探测定理”里的(x),考虑(a^{frac{p-1}{2}})(mod P)意义下的值:

  • 如果既不是(1)也不是(P-1),说明(P)不是质数。直接反回false

  • 如果是(P-1),那么不能利用二次探测继续递归,说明目前无法验证(P)为合数,返回true

  • 如果是(1),则我们可以递归下去,令(p:=frac{p-1}{2}+1)。继续判断(a^{frac{p-1}{2}})

如果每一层递归都做一次快速幂,复杂度会多个(log)。事实上,我们可以先把(p-1)每次除以(2)直到不能除为止,这相当于递归的最底层。然后在最底层做一次快速幂。再从底向上逐层推出当前层要判定的数的值。每一层的值都是在它下一层的基础上平方一下。求出每一层的数后,把它们存在一个vector里,从后往前扫一遍,就能模拟出递归的过程了。

复杂度(O(log p))

参考代码:

inline ll mul_mod(ll a,ll b,ll m=MOD) {
	a%=m;b%=m;
	ll c=(long double)a*b/m;
	ll ans=a*b-c*m;
	if(ans<0)ans+=m;
	if(ans>=m)ans-=m;
	return ans;
}
inline ll pow_mod(ll x,ll i,ll m=MOD){
	ll y=1;
	while(i){
		if(i&1LL)y=mul_mod(y,x,m);
		x=mul_mod(x,x,m);
		i>>=1;
	}
	return y;
}

const int B[5]={2,3,5,7,11};
bool check(ll x,int b){
	if(!(x&1))return false;
	ll k=x-1;
	while(!(k&1LL))k>>=1;
	ll t=pow_mod(b,k,x);
	vector<ll>vec;
	while(k!=x-1)k<<=1,t=mul_mod(t,t,x),vec.pb(t);
	if(t!=1)return false;
	for(int i=vec.size()-2;i>=0;--i){
		if(vec[i]!=1&&vec[i]!=x-1)return false;
		if(vec[i]==x-1)return true;
	}
	return true;
}
bool is_prime(ll x,int RC=2){
	/*
	RC 容错
	是一个1~5的整数,默认为2
	数字越大,准确度越高,效率越低
	*/
	if(x==1)return false;
	for(int i=0;i<RC;++i)if(x==B[i])return true;
	for(int i=0;i<RC;++i)if(!check(x,B[i]))return false;
	return true;
}

解决了素数判定,CASE 8~10的代码就相当简单了:

void solve5(){
	int n=read();while(n--){
		ll l=readll(),r=readll();
		for(ll i=l;i<=r;++i){
			if(is_prime(i,4))putchar('p');
			else putchar('.');
		}
		putchar('
');
	}
}

CASE (11sim 13)

套路和CASE 8~10差不多,只不过之前是对区间每个数做素数判定,现在是求区间每个数的莫比乌斯函数

(x=prod_{i=1}^{m}p_i^{c_i}),则莫比乌斯函数(mu(x))定义为:

[mu(x)=left{egin{matrix} 0&&{exists c_igeq 2}\ (-1)^m&& ext{otherwise}\ end{matrix} ight. ]

与测试点(8)类似,测试点(11)我们可以用线性筛法轻松解决。

对于测试点(12,13)(l,rleq 10^{18}),需要一点数论技巧。

我们先求出([1,10^6])的所有质数,用它们筛掉([l,r])的数的所有(leq 10^6)的质因子,筛的过程中顺便维护一下每个数的(mu)值:每新增一个质因子,(mu)值要乘以(-1),若一个质因子出现了大于(1)次,则(mu)值直接等于(0)

每个数被筛去这些小质因子后,留下的部分的质因子数最多不超过2。具体来讲,分如下几种情况:

  • 没有其他质因子,即留下的数(=1)

  • 留下的数是一个质数。此时(mu)值要乘以(-1)

  • 留下的数是一个完全平方数,那它肯定是某个质数的平方。此时(mu)的值直接(=0)

  • 留下的数是两个不同的质数的积,(mu)值不变。

枚举每个(leq 10^6)的质因子,再枚举其倍数,这样做的复杂度是调和级数,即(O(nlog n))级别。输出答案的部分,对每个数可能要使用大素数判定,复杂度(O(nlog p))(p)是值的大小,在(10^{18})级别。常数很大,我朴素地实现了一下,在本地要跑6秒。做了如下优化后勉强跑到2秒左右:

  • 筛质数时筛到某个(mu)值已经为(0)的数,直接跳过。

  • 最后的输出答案部分中,如果(mu)值已经为(0),则不要再进行素数判定。

参考代码:

using Math::is_prime;
using Math::sieve_mu;
ll val[1000006];
int mu[1000006];
void solve6(){
	#define C(x) ((x==0)?('0'):((x==1)?'+':'-'))
	sieve_mu();
	int n=read();
	//cout<<Math::cnt<<" "<<Math::p[78499]<<endl;
	for(int i=1;i<=n;++i){
		ll l=readll(),r=readll();
		if(l<=1000000){
			for(ll j=l;j<=min(r,1000000LL);++j)putchar(C(Math::mu[j]));
		}
		if(r<=1000000){
			putchar('
');
			continue;
		}
		l=max(l,1000001LL);
		for(ll j=l;j<=r;++j)val[j-l]=j,mu[j-l]=1;
		for(int j=1;j<=Math::cnt;++j){
			ll x=l/Math::p[j],y=r/Math::p[j];
			for(ll k=x;k<=y;++k){
				if(k*Math::p[j]<l)continue;
				//assert(k*Math::p[j]>=l);
				//assert(k*Math::p[j]<=r);
				if(!mu[k*Math::p[j]-l])continue;
				int cnt=0;
				while(val[k*Math::p[j]-l]%Math::p[j]==0){
					val[k*Math::p[j]-l]/=(ll)Math::p[j];
					cnt++;
					if(cnt>1){mu[k*Math::p[j]-l]=0;break;}
				}
				if(cnt==1)mu[k*Math::p[j]-l]*=-1;
			}
		}
		for(ll j=l;j<=r;++j){
			//cout<<val[j-l]<<" "<<mu[j-l]<<endl;
			if(val[j-l]==1||mu[j-l]==0)putchar(C(mu[j-l]));
			else if(is_prime(val[j-l]))putchar(C(-mu[j-l]));
			else{
				ll s=sqrt(val[j-l]);
				if(s*s==val[j-l])putchar(C(0));
				else putchar(C(mu[j-l]));
			}
		}
		putchar('
');
	}
	#undef C
}

CASE (14sim 16)

根据前面的经验,结合功能编号g,不难想到这部分数据是要让大家判断区间([l,r])里每个数是不是(m)原根


简单讲一下原根的定义和性质。

根据欧拉定理我们知道,(a,m)互质时,(a^{varphi(m)}equiv 1pmod m)。所以显然存在一类(x)使得(a^xequiv 1pmod m)(如(x=varphi(m)))。但(varphi(m))不一定是最小的(x)。我们希望找到满足条件的最小的正整数(x),即(min{x|x>0,a^xequiv 1pmod m}),并将这样的(x)命名为(a)(m)的阶。

对于一个(a),若(a)(m)的阶为(varphi(m)),则称(a)(m)的原根。

结论1:设(a)(m)的阶为(y),那么({x|x>0,a^xequiv 1pmod m})这个集合一定等于({ky|kinmathbb{N_+}}),也即,(a^xequiv 1pmod mLeftrightarrow y|x)

证明:

首先所有(y)的倍数(x)一定满足(a^xequiv 1pmod m)。因为根据定义,在(mod m)意义下,每乘以(y)(a)就会形成一个周期。

考虑如果存在一个满足(a^xequiv 1pmod m)(x)不被(y)整除,设(x=qy+r(0<r<y))。因为(a^{qy}equiv 1pmod m),所以(a^requiv 1pmod m),这与(y)(a)(m)的阶矛盾。故原命题不成立。即:所有满足(a^xequiv 1pmod m)(x)都能被(y)整除。

结论2:一个数(m)有原根,当且仅当(m)可以表示为如下四种形式的任意一种:(2,4,p^k,2p^k),其中(p)是一个质数。

证明略。

结论3:如果(m)有一个原根(g),那么(g,g^2,dots,g^{varphi(m)})(mod m)意义下的值一定恰好通过(m)的简化剩余系。

证明:

首先我们要证明,对任意(1leq j<kleqvarphi(m))(g^j eq g^kpmod m)。反证法。假设(g^kequiv g^jpmod m),则(frac{g^k}{g^j}=g^{k-j}equiv1 pmod m)。即:存在一个(1<k-j<varphi(m))使得(g^{k-j}equiv 1pmod m),这与(g)(m)的原根矛盾。

于是我们证明了我们取的数在(mod m)意义下两两不同。又因为我们取了(varphi(m))个数,所以这些数恰好通过(m)的简化剩余系。


回到本题,我们要判断一个数(a)是否是(m)的原根。朴素的算法我们要枚举(a)(1simvarphi(m))次幂,复杂度是(O(varphi(m)))的,无法承受。

但是根据结论1,任何一个满足(a^xequiv 1pmod m)(x)都是(a)(m)的阶的倍数。(varphi(m))显然是其中一个满足条件的(x)。要找(a)(m)的阶,我们只需要枚举(varphi(m))的约数即可。

更进一步,如果一个约数(x)满足(a^xequiv 1pmod m),那么(x)的倍数也满足。于是我们只需要枚举(varphi(m))的每个质因数(p_i),看(frac{varphi(m)}{p_i})是否满足条件。如果存在(frac{varphi(m)}{p_i})满足条件,说明(a)(m)的阶一定小于(varphi(m))(a)必不是(m)的原根。否则,可以说明(a)必是(m)的原根。


对于测试点(14),可以预处理出(varphi(m)=998244352)的所有质因数,判断每个(ain[l,r])是否是原根时都扫一遍所有质因数。复杂度(O(sqrt{m}+ntlog m)),其中(n=r-l+1)(t)(varphi(m))的质因数个数,(log m)是判断时做快速幂的复杂度。

测试点(15)的区间长度很大,但模数(m)较小,考虑换一种思路。先求出(m)的一个原根(g)。因为(m)是质数,那么根据结论3([1,m-1])的每个数都能被表示为(g^i(1leq i<m))。设当前要判断的数(a=g^x)。若(gcd(x,m-1) eq 1),则(a)不是原根。因为如果(gcd(x,m-1)=d),则(a^{frac{m-1}{d}}=g^{frac{x(m-1)}{d}}equiv 1pmod m)。否则不存在(m-1)的约数(d),使得(a^{frac{m-1}{d}}equiv 1pmod m),根据前面基于结论1的讨论,此时(a)一定是原根。

如果每次暴力求(gcd)显然还是太慢了,我们可以用(m-1)的所有质因数依次筛去它们的所有倍数。具体实现见代码。

测试点(16)里的?,我们从(10^9)开始暴力求。注意,根据结论2以及基于对出题人的信任,此处我们还是只枚举质数就好。用下面的代码大约10分钟左右可以求出来。

void get_pr(){
	int L=1e9,R=2e9,l=233333333,r=234133333;
	freopen("software16.ans","r",stdin);
	string str;cin>>str;cin>>str;cin>>str;
	for(int P=L,c=0;P<=R;++P){
		if(!is_prime(P))continue;
		++c;
		if(c%10000==0)cerr<<P<<endl;
		div(P-1);
		bool fail=0;
		for(int i=l;i<=r;++i){
			char ch=(is_pr(i,P)?'g':'.');
			if(ch!=str[i-l]){
				fail=1;break;
			}
		}
		if(fail)continue;
		cerr<<"! "<<P<<endl;
		return;
	}
}

CASE 14~16的代码:

int p[100005],cnt;
void div(int x){
    cnt=0;
    for(int i=2;i*i<=x;++i)if(x%i==0){
        while(x%i==0)x/=i;
        p[++cnt]=i;
    }
    if(x!=1)p[++cnt]=x;
}
inline int is_pr(int x,int P){
    for(int i=1;i<=cnt;++i){
        if(pow_mod(x,(P-1)/p[i],P)==1)return false;
    }
    return true;
}
bool vis[20000007];
int ind[20000007];
void solve7(){
	div(998244353-1);
	int n=read();while(n--){
		int l=read(),r=read(),P;
		if(l==233333333)P=1515343657,div(P-1);
		else P=read();
		if(P==998244353||P==1515343657){
			for(int i=l;i<=r;++i){
				if(is_pr(i,P))putchar('g');
				else putchar('.');
			}
		}
		else{
			div(P-1);
			for(int i=1;i<=cnt;++i){
				for(int j=1;j*p[i]<=P-1;++j){
					vis[j*p[i]]=1;
				}
			}
			int g=0;
			for(int i=1;i<P;++i)if(is_pr(i,P)){g=i;break;}
			for(int i=1,t=g;i<=P-1;++i,t=(ll)t*g%P)ind[t]=i;//对数
			for(int i=l;i<=r;++i){
				if(vis[ind[i]])putchar('.');
				else putchar('g');
			}
		}
		putchar('
');
	}
}

完整代码

原文地址:https://www.cnblogs.com/dysyn1314/p/12371784.html