题解 LOJ3050 「十二省联考 2019」骗分过样例
CASE \(1\sim 3\)
\(n\)组测试数据，每次输入一个数\(x\)，求\(19^x\)。
测试点\(1\)，\(x=0,1,\dots n-1\)，可以直接递推。
测试点\(2\)要开long long并用快速幂。
测试点\(3\)中\(x\)超出了long long范围。根据欧拉定理，当\(a,p\)互质时，\(a^b\equiv a^{b\bmod\varphi(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\n",pow_mod(19,x));
}
}
CASE \(4\sim 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_x\cdot19^{y-x}\equiv ans_y\pmod p\)。
所以模数\(p\)就一定是\(ans_x\cdot 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 \(6\sim 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\n",ans[x]);
else printf("%d\n",ans[(x-ST)%LEN+ST]);
}
}
CASE \(8\sim 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 1\pmod p \]
我们选取一些\(a\)，然后把\(p\)带入，进行判定。若对于其中某个\(a_i\)，\(p\)不满足费马小定理，则可以确定\(p\)不是质数。
但这个条件是必要不充分条件。也就是说，存在一类伪素数，它们虽然不是素数，但满足费马小定理。
为避免这种情况，我们的算法要进行二次探测。它基于如下定理：
对一个素数\(p\)，若有\(x^2\equiv 1\pmod p\)，则\(x\equiv 1\pmod p\)或\(x\equiv p-1\pmod p\)。
证明：移项得\(x^2-1\equiv 0\pmod p\)。用平方差公式展开得\((x+1)(x-1)\equiv 0\pmod p\)。所以\(x\equiv 1\pmod p\)或\(x\equiv p-1\pmod p\)。
很遗憾，这个定理同样是必要不充分条件。但我们可以利用它来加强我们的条件，提高素数判定的准确性。
请注意，接下来我们讨论的是一个递归的过程，在递归过程中，指数\(p\)是在不断变化的，而模数\(P\)是始终不变的。初始时，\(P=p\)，都是我们要判定的那个数。
首先，经过费马小定理的验证，我们已经确定\(a^{p-1}\equiv 1\pmod P\)。如果\(p-1\)是偶数，则有\(\left(a^{\frac{p-1}{2}}\right)^2\equiv 1\pmod P\)。那么可以将\(a^{\frac{p-1}{2}}\)视作“二次探测定理”里的\(x\)，考虑\(a^{\frac{p-1}{2}}\)在\(\bmod 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('\n');
}
}
CASE \(11\sim 13\)
套路和CASE 8~10差不多，只不过之前是对区间每个数做素数判定，现在是求区间每个数的莫比乌斯函数。
设\(x=\prod_{i=1}^{m}p_i^{c_i}\)，则莫比乌斯函数\(\mu(x)\)定义为：
与测试点\(8\)类似，测试点\(11\)我们可以用线性筛法轻松解决。
对于测试点\(12,13\)，\(l,r\leq 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(n\log n)\)级别。输出答案的部分，对每个数可能要使用大素数判定，复杂度\(O(n\log 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('\n');
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('\n');
}
#undef C
}
CASE \(14\sim 16\)
根据前面的经验，结合功能编号g，不难想到这部分数据是要让大家判断区间\([l,r]\)里每个数是不是\(m\)的原根。
简单讲一下原根的定义和性质。
根据欧拉定理我们知道，\(a,m\)互质时，\(a^{\varphi(m)}\equiv 1\pmod m\)。所以显然存在一类\(x\)使得\(a^x\equiv 1\pmod m\)（如\(x=\varphi(m)\)）。但\(\varphi(m)\)不一定是最小的\(x\)。我们希望找到满足条件的最小的正整数\(x\)，即\(\min\{x|x>0,a^x\equiv 1\pmod m\}\)，并将这样的\(x\)命名为\(a\)对\(m\)的阶。
对于一个\(a\)，若\(a\)对\(m\)的阶为\(\varphi(m)\)，则称\(a\)是\(m\)的原根。
结论1：设\(a\)对\(m\)的阶为\(y\)，那么\(\{x|x>0,a^x\equiv 1\pmod m\}\)这个集合一定等于\(\{ky|k\in\mathbb{N_+}\}\)，也即，\(a^x\equiv 1\pmod m\Leftrightarrow y|x\)。
证明：
首先所有\(y\)的倍数\(x\)一定满足\(a^x\equiv 1\pmod m\)。因为根据定义，在\(\bmod m\)意义下，每乘以\(y\)个\(a\)就会形成一个周期。
考虑如果存在一个满足\(a^x\equiv 1\pmod m\)的\(x\)不被\(y\)整除，设\(x=qy+r(0<r<y)\)。因为\(a^{qy}\equiv 1\pmod m\)，所以\(a^r\equiv 1\pmod m\)，这与\(y\)是\(a\)对\(m\)的阶矛盾。故原命题不成立。即：所有满足\(a^x\equiv 1\pmod m\)的\(x\)都能被\(y\)整除。
结论2：一个数\(m\)有原根，当且仅当\(m\)可以表示为如下四种形式的任意一种：\(2,4,p^k,2p^k\)，其中\(p\)是一个质数且\(p\geq 3\)。
证明略。
结论3：如果\(m\)有一个原根\(g\)，那么\(g,g^2,\dots,g^{\varphi(m)}\)在\(\bmod m\)意义下的值一定恰好通过\(m\)的简化剩余系。
证明：
首先我们要证明，对任意\(1\leq j<k\leq\varphi(m)\)，\(g^j\neq g^k\pmod m\)。反证法。假设\(g^k\equiv g^j\pmod m\)，则\(\frac{g^k}{g^j}=g^{k-j}\equiv1 \pmod m\)。即：存在一个\(1<k-j<\varphi(m)\)使得\(g^{k-j}\equiv 1\pmod m\)，这与\(g\)是\(m\)的原根矛盾。
于是我们证明了我们取的数在\(\bmod m\)意义下两两不同。又因为我们取了\(\varphi(m)\)个数，所以这些数恰好通过\(m\)的简化剩余系。
回到本题，我们要判断一个数\(a\)是否是\(m\)的原根。朴素的算法我们要枚举\(a\)的\(1\sim\varphi(m)\)次幂，复杂度是\(O(\varphi(m))\)的，无法承受。
但是根据结论1，任何一个满足\(a^x\equiv 1\pmod m\)的\(x\)都是\(a\)对\(m\)的阶的倍数。\(\varphi(m)\)显然是其中一个满足条件的\(x\)。要找\(a\)对\(m\)的阶，我们只需要枚举\(\varphi(m)\)的约数即可。
更进一步，如果一个约数\(x\)满足\(a^x\equiv 1\pmod 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\)的所有质因数，判断每个\(a\in[l,r]\)是否是原根时都扫一遍所有质因数。复杂度\(O(\sqrt{m}+nt\log m)\)，其中\(n=r-l+1\)，\(t\)是\(\varphi(m)\)的质因数个数，\(\log m\)是判断时做快速幂的复杂度。
测试点\(15\)的区间长度很大，但模数\(m\)较小，考虑换一种思路。先求出\(m\)的一个原根\(g\)。因为\(m\)是质数，那么根据结论3，\([1,m-1]\)的每个数都能被表示为\(g^i(1\leq i<m)\)。设当前要判断的数\(a=g^x\)。若\(gcd(x,m-1)\neq 1\)，则\(a\)不是原根。因为如果\(gcd(x,m-1)=d\)，则\(a^{\frac{m-1}{d}}=g^{\frac{x(m-1)}{d}}\equiv 1\pmod m\)。否则不存在\(m-1\)的约数\(d\)，使得\(a^{\frac{m-1}{d}}\equiv 1\pmod 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('\n');
}
}
