天下一プログラマーコンテスト2015 予選B : E - 天下一演算

これは自力では思いつかないわ…。
http://tenka1-2015-qualb.contest.atcoder.jp/tasks/tenka1_2015_qualB_e

問題

M未満の非負整数で構成されるN要素の数列Aがある。
この数列に対し、全要素を10倍する処理を何度か行った後、任意の要素を1加算する処理を何度か行った結果、全要素がMの倍数になった。

条件を満たすような、全要素10倍と1加算の処理回数の和の最小値を求めよ。

解法

処理の過程において、各要素のMの剰余を考える。

各要素に10倍を繰り返す場合、10とMが互いに素なら、φ(M)の約数で要素は周期的に変化する。

10とMが互いに素でない場合でも、16回以上10倍すれば、以後φ(M)の約数で要素は周期的に変化する。
Mは10^5以下なので、10倍を16回行うと、その要素(のMの剰余)とMは、素因数分解したときの2と5の指数は一致する(2^17＞10^5のため)。
A[i]を10^16した値をBとし、G=GCD(10^16,M)とすると、B/GとM/Gは互いに素なので、その周期はφ(M/G)となるためである。

よって、まず最初に10倍を16回まで行うケースにおける処理回数は愚直に総当たりして求めてしまおう。
あとはA[i]=A[i]*10^16%MでA[i]を置き換えれば、各A[i]は周期的であるとみなせる。

まず1～(M-1)の各値について、A中に何度登場するかを求め、xが登場する回数をB(x)とする。
(A[i]=0の場合、10倍を何回してもそれ以上1加算は不要なので無視してよい)

あるi=1～(M-1)において、C=[i,(i*10^1)%M,(i*10^2)%M, ..., (i*10^(P-1))%M]という数列を考える。(Pは周期、すなわち(i*10^P)%M==iとする)
という多項式と、という多項式を考える。
するとにおける次の項の係数は、j回10倍処理を行ったとき、A[i]のうちCに含まれる要素における10倍処理後の1加算回数に一致する。
よってFFTを用いてを求めることで、10倍処理回数に応じた1加算回数を求めることができる。

Cは1～(M-1)を全てを含むとは限らないので、1～(M-1)を全部1回ずつ含むよう、いくつかのCにおいて上記FFTを行い、それぞれ10倍処理後の1加算回数を求めることで、A[i]全体を10倍処理後における1加算回数を求めることができる。

Mの値によっては、Cが短周期の複数の数列となりうるため、FFTを行う場合は数列数が少ない場合は計算対象の要素数を減らす等、計算量を落とす工夫が必要。

int N,M;
int A[101010];
int num[101010];
int used[101010];
ll tret[201010];
ll ret=1LL<<60;
typedef complex<double> Comp;

vector<Comp> fft(vector<Comp> v, bool rev=false) {
	int n=v.size(),i,j,m;
	
	for(i=0,j=1;j<n-1;j++) {
		for(int k=n>>1;k>(i^=k);k>>=1);
		if(i>j) swap(v[i],v[j]);
	}
	for(int m=2; m<=n; m*=2) {
		double deg=(rev?-1:1) * 2*acos(-1)/m;
		Comp wr(cos(deg),sin(deg));
		for(i=0;i<n;i+=m) {
			Comp w(1,0);
			for(int j1=i,j2=i+m/2;j2<i+m;j1++,j2++) {
				Comp t1=v[j1],t2=w*v[j2];
				v[j1]=t1+t2, v[j2]=t1-t2;
				w*=wr;
			}
		}
	}
	if(rev) FOR(i,n) v[i]*=1.0/n;
	return v;
}

vector<Comp> MultPoly(vector<Comp> P,vector<Comp> Q,bool resize=false) {
	if(resize) {
		int s=2;
		while(s*2<P.size()+Q.size()) s*=2;
		P.resize(s*2);Q.resize(s*2);
	}
	P=fft(P), Q=fft(Q);
	for(int i=0;i<P.size();i++) P[i]*=Q[i];
	return fft(P,true);
}

void solve() {
	int i,j,k,l,r,x,y; string s;
	
	cin>>M>>N;
	FOR(i,N) cin>>A[i];
	FOR(i,20) {
		ll tot=i;
		FOR(x,N) tot+=(M-A[x])%M, A[x]=A[x]*10%M;
		ret=min(ret,tot);
	}
	
	int alpha=__gcd(1<<20,M)*__gcd(125*125*125,M);
	if(alpha==M) return _P("%lld\n",ret);
	int per;
	
	FOR(i,N) num[A[i]]++;
	
	vector<vector<ll>> loops;
	for(i=alpha;i<M;i+=alpha) if(used[i]==0) {
		vector<Comp> nu,ad;
		int per;
		x = i;
		FOR(per,M) {
			used[x]=1;
			ad.push_back(M-x);
			nu.push_back(num[x]);
			x=x*10%M;
			if(x==i) break;
		}
		per++;
		reverse(ALL(nu));
		
		FOR(j,per) ad.push_back(ad[j]);
		auto prod = MultPoly(ad,nu,true);
		
		vector<ll> lp;
		FOR(j,per) lp.push_back((ll)(prod[per-1+j].real()+0.1));
		loops.push_back(lp);
	}
	
	int maxper=1;
	FORR(r,loops) maxper=maxper/__gcd(maxper,(int)r.size())*r.size();
	FOR(i,maxper) {
		ll tot=20+i;
		FORR(r,loops) tot += r[i%r.size()];
		ret=min(ret,tot);
	}
	
	cout<<ret<<endl;
	
}

まとめ

よくこんなの本番中に思いつくな…。