TopCoder SRM 657 Div1 Medium PolynomialGCD

1ミスしたものの、自力で回答できた。
http://community.topcoder.com/stat?c=problem_statement&pm=13774

問題

はxがどんな整数でも6の倍数になることが知られている。
一般化してN要素の整数列S[i]が与えられたとき、という形の関数を考える。

xがどんな整数であっても、P(x)がyの倍数になるというような最大のyを求めよ。

解法

N以下の各素数pに対し、最小でpの何乗が含まれるかを考える。
たとえばN==10の時、p=3としてP(x)は以下の３値のうち最小値をqとするとyはp^qの倍数になる。

(注意：以下の数値は厳密には異なる)

• xが3の倍数の時、S[0]+S[3]+S[6]+S[9]
• xが(3の倍数+1)の時、S[1]+S[4]+S[7]
• xが(3の倍数+2)の時、S[2]+S[5]+S[8]

注意点として、xが3の倍数の時、x,(x-3),(x-6),(x-9)のどれかは9の倍数でもあり得る。
よって実際は先頭のケースではさらに以下のようにp^2を考慮し、以下の最小値を求めなければならない。

• xが9の倍数の時 2*S[0]+S[3]+S[6]+2*S[9]
• xが(9の倍数+3)の時 S[0]+2*S[3]+S[6]+S[9]
• xが(9の倍数+6)の時 S[0]+S[3]+2*S[6]+S[9]

このようにp^rの倍数となる項がp個以上ある場合、そのうち1/p個はp^(r+1)の倍数にもなることに注意してpが何回掛算されるか最小値を数えあげていく。
以下のコードでは、iがp^rの倍数である場合を総当たりし、周囲のp^t(t＜r)の倍数を数え上げている。

ll mo=1000000007;
const int prime_max = 10000;
int NP,prime[100000],divp[prime_max];

void cprime() {
	NP=0;
	ZERO(divp);
	for(int i=2;i<prime_max;i++) if(divp[i]==0) {
		prime[NP++]=i;
		for(int j=i*i;j>i&&j<prime_max;j+=i) divp[j]=i;
	}
}

class PolynomialGCD {
	public:
	int val(int p,string& s) {
		int i,N=s.size();
		int mi=1000000000;
		
		if(p>N) return 0;
		FOR(i,N) {
			int tot=0;
			int cand=100000;
			for(int P=p;cand>=p;P*=p) {
				cand=1;
				tot+=s[i];
				for(int d=P;d<N;d+=P) {
					if(i+d<N) tot += s[i+d], cand++;
					if(i-d>=0) tot += s[i-d], cand++;
				}
			}
			mi=min(mi,tot);
		}
		return mi;
	}
	
	int gcd(string s) {
		int i,x;
		cprime();
		FORR(r,s) r-='0';
		
		ll ret=1;
		FOR(i,NP) {
			int num=val(prime[i],s);
			FOR(x,num) ret = ret*prime[i]%mo;
		}
		return ret;
	}
}

まとめ

問題の発想が面白かった。