Digit (AOJ2392)

AOJ-ICPC

解きました 以下ネタバレ

L進法でL-1のことを便宜的に9と書きます
L(0)=1,L(1)=10,L(2)=19,L(3)=199,L(4)=199..99(9が2L+1個)...
みたいな感じに19...9になるのが簡単に帰納法でわかります
L(N)の9の個数b_Nの漸化式はb_N=1+L+L^2+...+L^(b_{N-1}-1)となる.
また、答えはb_{N+1}*(L-1)+1 % Mになる.

で,これを求めるのに,単純にb_i mod Mとかを求めても当然ダメで,それぞれのiに対してb_i mod M_iを求めて(M_{N+1}=問題文のM)順番に計算していく,ということをする必要がある.以下M_{i+1}からM_{i}を求めるとして,M_{i+1}のことをM,M_{i}のことをmと書く.

漸化式の形から,1+L+...L^(x-1)=0(mod M)となれば,b_{i}をmod xで求めれば良くなる(本当はLとM_iが互いに素でない時単純にこれではダメだけど後述)
このxの必要条件として,両辺に1+LをかければL^x=1(mod M)となる.とりあえずこれを満たす最小のxをoと置く.oはbaby-giantで求める.この時1+L+..+L^(o-1)=y(mod M)とおくと,求めるx(つまりm)はo*M/gcd(y,M)となる(xがoの倍数なのは必要条件で,かつoのところから足されるものはループを起こすので).
M_{i+1}からM_iを計算する時に,どんどん大きくなっていくのでは?とか思うかもしれないが,大きくならないことが示せる(数論)
で、ここは怪しいんだけど,M_{i+1}=M_iになる(その後ずっとM_iのままになる)か,そうでない時は結構M_iは小さくなるんじゃないか?という感じがしたのでそれを仮定したら通ってしまった(計算量解析してなし),実際どのテストケースでも値が異なるM_iは10種類程度に抑えられていた.
LとMが互いに素じゃない時は,始めの数項(具体的にはgcd(L^a,M)とgcd(L^(a+1),M)が同じになる時1+L+..+L^(a-1)まで(例えば,L=10,M=24のとき,1+10+4))を計算しておく,そっから先はループが出来る(上のgcdの値をgとおくと,gで合同式を割るとM'(=M/g)とLが互いに素になるので)ので,互いに素な時に帰着できた.ここらへんの細かい計算が面倒なだけです(あと,"始めの数項"の途中で計算が途切れてる場合(これはb_{i+1}が小さいと起こる)があるけど,b_{i}を適当に計算すればわかるけどb_4くらいまで手元でmodを取らずに求めとけば大丈夫です)
完全に説明がめんどくなった.以下コード.

#include <bits/stdc++.h>
#define rep(i,n) for(int i=0;i<(int)(n);i++)
#define rep1(i,n) for(int i=1;i<=(int)(n);i++)
#define all(c) c.begin(),c.end()
#define pb push_back
#define fs first
#define sc second
#define show(x) cout << #x << " = " << x << endl
#define chmin(x,y) x=min(x,y)
#define chmax(x,y) x=max(x,y)
using namespace std;
typedef long long ll;
ll a[100010];
ll b[100010];
ll c[100010];
ll s[100010];
ll m[100010];
ll gcd(ll x,ll y){
	if(y==0) return x;
	return gcd(y,x%y);
}
ll ex(ll x,ll p,ll mod){
	ll a=1;
	while(p){
		if(p%2) a=a*x%mod;
		x=x*x%mod;
		p/=2;
	}
	return a;
}
ll geo(ll x,ll p,ll mod){	//1+x+..+x^(p-1)
	if(p==0) return 0;
	ll a=0,s=1,t=1;
	while(p){
		if(p%2){
			a=(a+s*t)%mod;
			s=s*x%mod;
		}
		t=t*(1+x)%mod;
		x=x*x%mod;
		p/=2;
	}
	return a;
}
typedef pair<int,int> P;
ll getord(ll M,ll a){
	ll H=sqrt(M)+1;
	vector<P> babies;
	ll p=1;
	rep(i,H){
		babies.pb(P(p,i));
		p=p*a%M;
	}
	ll q=p;
	sort(all(babies));
	rep1(i,H){
		int pos=upper_bound(all(babies),P(p,H))-babies.begin()-1;
		if(pos>=0&&babies[pos].fs==p) return i*H-babies[pos].sc;
		p=p*q%M;
	}
}
ll pre,parg;
void getnext(ll M,ll L,int id){
	ll p=L;
	ll gs[101];
	gs[0]=1;
	rep(i,100){
		gs[i+1]=gcd(p,M);
		if(gs[i]==gs[i+1]){
			a[id]=i;
			s[id]=geo(L,i,M);
			ll g=gs[i];
			c[id]=ex(L,i,M);
			M/=g;
			break;
		}
		p=p*L%M;
	}
	if(M==pre){
		m[id]=parg;
		return;
	}
	ll o=getord(M,L);
	ll ss=geo(L,o,M);
	m[id]=M/gcd(M,ss)*o;
	pre=M,parg=m[id];
}
ll getans(ll b,ll M,ll L){
	return (2*ex(L,b,M)+M-1)%M;
}
int solve(ll N,ll M,ll L){
	pre=-1;
	if(N==0) return 1;
	if(N==1) return L%M;
	if(N==2) return getans(1,M,L);
	if(N==3) return getans(2,M,L);
	if(N==4) return getans(2*(L+1),M,L);
	N++;
	m[N]=M;
	for(ll x=N-1;x>=4;x--){
		getnext(m[x+1],L,x);
	}
	b[4]=2*(L+1)%m[4];
	for(ll x=4;x<N;x++){
		ll p=(b[x]-a[x])%m[x];
		if(p<0) p+=m[x];
		b[x+1]=2*(s[x]+c[x]*geo(L,p,m[x+1]))%m[x+1];
	}
	return (b[N]*(L-1)+1)%M;
}
int main(){
	for(int t=1;;t++){
		ll N,M,L;
		cin>>N>>M>>L;
		if(M==0) break;
		printf("Case %d: %d\n",t,solve(N,M,L));
	}
}