CodeChef-MISINT2 : Misinterpretation 2

codechef

問題概要

長さnの文字列で、偶数番目+奇数番目と連結した文字列ともとの文字列が一致するのは何通りあるかの和を区間[L,R]について求めよ。L-R<5*10^4, R<10^10。

解法

小さいケースについてはunion-findでつなげて親の個数を数えてやればよい。これで試すと偶数項と奇数項の間に分かりやすい関係が見つかり片方を数列辞典に投げると見事にヒットする。ord_2(n)を2^mが1とmod nで等しくなるような最小の正整数mで定めると、奇数項についてa[n]=Σ_{d|n}phi(d)/ord_2(d)となるらしい。制約が明らかに区間篩を臭わせていて、実際に区間篩を行ってついでに素因数分解しておくことによってphi(d)はdを列挙するのと同時に計算できる。ord_2(Π p_i ^ e_i) = LCM(ord_2(p_i^e_i))が成り立つこと予想できるのでこれを使ってord_2(p_i)もそれなりの速度で計算できる。小さいところでord_2(p^k) = p * ord_2(p^(k-1))がk>=2について成り立つらしいことが確かめられる(証明分からない)のでそれを使って高速化したりする。k=1のときは、phi(p)=p-1を素因数分解(ord_2(n) | phi(n)なので)して愚直に求めて小さいところでメモ化したりしてやると良い。
この問題、ord_2の計算にphi(d)の約数を調べるだけでも頑張って高速化すれば想定解の2倍程度の速度までは平気で出るので、必死に高速化頑張ったけど結果的にはそんなに意味はなかった。高速化の練習にはなったと思いたい。あと上限10^10が微妙にいやらしくてオーバーフローに注意する必要もある。

acceptされたコード

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

typedef unsigned long long int64;

const int MAX_RT = 1e5+10;
const int MAX_W = 5e4+10;
const int64 MOD = 1e9 + 7;
const int MEMO_SIZE = 1e5;
const int MAX_SIZE = 1e5+100;
const int MAX_PRIME = 1e4;


int64 R, L;
int fs[MAX_W];
int64 ps[MAX_W][11];
int ns[MAX_W][11];
int64 pps[11];
int es[11];

bool visited[MAX_RT][30];
int64 memo[MAX_RT][30];

bool vis2[MEMO_SIZE];
int64 memo2[MEMO_SIZE];

int64 bins[40];

bool isPrime[MAX_SIZE];
int P;
int primes[MAX_PRIME];
int64 step[MAX_PRIME];

int phiF;
int64 phipps[11];
int phinns[11];

int64 prev_n, mem;

void seive(){
	memset(isPrime, true, sizeof(isPrime));
	isPrime[0] = isPrime[1] = false;
	for(int i=2; i*i<MAX_SIZE; i++)if(isPrime[i]){
		for(int j=i<<1; j<MAX_SIZE; j+=i){
			isPrime[j] = false;
		}
	}
	for(int i=2, st=1e5+10; i<st; i++){
		if(isPrime[i]){
			step[P] = (int64)i*i;
			primes[P++] = i;
		}
	}
}

void init(){
	scanf("%lld%lld", &L, &R);
	memset(fs, 0, sizeof(fs));
}

int64 lcm(int64 x, int64 y){
	return x / __gcd(x, y) * y;
}

int64 pow_binary2(int64 n, int64 mod){
	int64 ret = 1;

	if(mod < (1LL<<32)){
		for(;n;n&=n-1){
			ret = ret * bins[__builtin_ctzll(n)] % mod;
		}
	}
	else{
		for(;n;n&=n-1){
			ret = (((ret * (bins[__builtin_ctzll(n)]>>16) % mod)<<16) + ret * (bins[__builtin_ctzll(n)]&0xffff)) % mod;
		}
	}

	return ret;
}

int64 pow_binary(int64 x, int64 n, int64 mod = MOD){
	int64 ret = 1;
	for(;n;n>>=1){
		if( n&1 ){
			ret = ret * x % mod;
		}
		x = x * x % mod;
	}

	return ret;
}

int64 pow_binary3(int64 x, int64 n, int64 mod){
	int64 ret = 1;

	if(mod < (1LL<<32)){
		for(;n;n>>=1){
			if( n&1 ){
				ret = ret * x % mod;
			}
			x = x * x % mod;
		}
	}
	else{
		for(;n;n>>=1){
			if( n&1 ){
				ret = (((ret * (x>>16) % mod)<<16) + ret * (x&0xffff))%mod;
			}
			x = (((x * (x>>16) % mod)<<16) + x * (x&0xffff))%mod;
		}
	}

	return ret;
}

int64 ord2pk(int64 p, int k){
	if(p < MAX_RT && visited[p][k]){
		return memo[p][k];
	}

	if(k == 1){
		int64 n = p;

		int64 m = p - 1;
		phipps[0] = 2;
		phinns[0] = __builtin_ctzll(m);
		phiF = 1;
		m >>= __builtin_ctzll(m);

		for(int i=1; step[i]<=m; i++){
			if(m%primes[i] == 0){
				phipps[phiF] = primes[i];
				phinns[phiF] = 0;
				for(;m%primes[i]==0;m/=primes[i]){
					phinns[phiF]++;
				}
				phiF++;
			}
		}
		if(m!=1){
			phipps[phiF] = m;
			phinns[phiF] = 1;
			phiF++;
		}

		bins[0] = 2;
		if(n > (1LL<<32)){
			for(int i=1, ed=64-__builtin_clzll(n); i<ed; i++){
				bins[i] = (((bins[i-1] * (bins[i-1]>>16) % n)<<16) + bins[i-1] * (bins[i-1]&0xffff))%n;
			}
		}
		else{
			for(int i=1, ed=64-__builtin_clzll(n); i<ed; i++){
				bins[i] = bins[i-1] * bins[i-1] % n;
			}
		}


		int64 phi = p - 1;
		int64 ret = phi;

		for(int i=0; i<phiF; i++){
			int64 tt = ret;
			for(int j=0; j<phinns[i]; j++){
				tt /= phipps[i];
			}

			int64 xx = pow_binary2(tt, n);
			if(xx == 1){
				ret = tt;
				continue;
			}
			for(int j=0; j<phinns[i]-1; j++){
				tt *= phipps[i];
				xx = pow_binary3(xx, phipps[i], n);
				if(xx == 1){
					ret = tt;
					break;
				}
			}
		}

		if(p < MAX_RT){
			visited[p][k] = true;
			return memo[p][k] = ret;
		}
		return ret;
	}

	else{
		int64 ret = ord2pk(p, 1);
		for(int i=0; i<k-1; i++){
			ret *= p;
		}

		if(p < MAX_RT){
			visited[p][k] = true;
			return memo[p][k] = ret;
		}
		return ret;
	}
}

int64 ord2(int F){

	int64 ret = 1;
	for(int i=0; i<F; i++)if(es[i] > 0){
		ret = lcm(ret, ord2pk(pps[i], es[i]));
	}

	return ret;
}

int64 ord(int64 n){
	int l = __builtin_ctzll(n+1);
	int64 m = (n+1)>>l, ret = l;
	for(;m!=1;){
		l = __builtin_ctzll(n+m);
		m = (n+m) >> l;
		ret += l;
	}

	return ret;
}

int64 phi_ord(int64 n, int64 phi, int depth){
	if(n >= MEMO_SIZE){
		return phi / ord2(depth);
	}

	if(vis2[n]){
		return memo2[n];
	}
	vis2[n] = true;
	return memo2[n] = phi / ord2(depth);
}

int64 dfs(int depth, int64 cur, int64 phi, int idx){
	if(depth == fs[idx]){
		return phi_ord(cur, phi, depth);
	}

	int64 ans = 0, p = ps[idx][depth], pp=p-1;
	pps[depth] = p;
	for(int i=0, ed = ns[idx][depth]; i<=ed; i++){
		es[depth] = i;
		if(i==0){
			ans += dfs(depth+1, cur, phi, idx);
		}
		else{
			cur *= p;
			ans += dfs(depth+1, cur, phi*pp, idx);
			pp *= p;
		}
	}
	return ans;
}

void check(int64 n){
	int idx = n - L, f = fs[idx];
	int64 m = n;
	for(int i=0; i<f; i++){
		int64 p = ps[idx][i];
		ns[idx][i] = 0;
		for(;m%p==0;m/=p){
			ns[idx][i]++;
		}
	}
	if(m!=1){
		ps[idx][f] = m;
		ns[idx][f] = 1;
		fs[idx]++;
	}
}

int64 calc(int64 n){
	if(prev_n == n){
		return mem;
	}

	if(n&1){
		check(n);
		prev_n = n;
		return mem = dfs(0, 1, 1, n-L);
	}
	return calc(n+1) - 1;
}

void prepare(){
	for(int i=0; i<P; i++){
		for(int64 j=max((int64)2, (L+primes[i]-1)/primes[i])*primes[i];j<=R+1; j+=primes[i]){
			ps[j-L][fs[j-L]] = primes[i];
			fs[j-L]++;
		}
	}

}

int solve(){
	int64 ans = 0;

	prepare();
	prev_n = 1LL<<60;
	for(int64 i=L; i<=R; i++){
		ans += pow_binary(26, calc(i));
		if(ans >= MOD){
			ans -= MOD;
		}
	}

	return (int)ans;
}

int main(){
	seive();

	int T;
	scanf("%d", &T);
	while(T--){
		init();
		printf("%d\n", solve());
	}

	return 0;
}