kmjp's blog

競技プログラミング参加記です

yukicoder : No.248 ミラー君の宿題

yukicoder

こんなの自分の知識では解説なしに解けません…。
http://yukicoder.me/problems/580

問題

N個の異なる奇素数p[i]が与えられる。
p[i]の積であるPを、(p[i]-1)の積であるQを用いて素因数分解したい。

f(X,Q)をXの素因数分解した結果を返す関数とする。
この関数は以下のように定義され、Pの約数を素因数分解するのに利用できる。

  1. X=1なら空集合を返す。
  2. Xが素数なら{X}を返す。
  3. 1～Xの範囲の整数を等確率で選び、その値をrとする。さらにgをGCD(r,X)とする。
    1. g==Xならf(X,Q)を再帰的に行う。
    2. 1＜g＜XならgはXの約数なので、f(g,Q)∪f(X/g,Q)を再帰的に行う。
    3. g==1の時、Q=2^e*q (qは奇数)と表現できるとして、c=r^p (mod X)を求め、c^2 != 1 (mod X)である限りc = c^2 (mod P)で更新する。
      1. c=1またはc=X-1ならf(X,Q)を再帰的に行う。
      2. 1＜c＜X-1ならhをGCD(c-1,X)とし、f(h,Q)∪f(X/h,Q)を再帰的に行う。

f(P,Q)を求める際、関数fを呼ぶ回数の期待値を求めよ。

解法

Writer解説そのままなので、概要はそちらを見た方が良いと思います。

Xがp[i]のうちどれの積であるか、を状態としてbitDPまたはメモ化再帰で処理する。

上記手順の3-2までは割と簡単。
というのもrに対し、1＜g＜Xとなる確率は簡単に求められる。
gはp[i]のうち一部a[i]の積であり、残りX/rはp[i]のうち一部b[i]の積であるとする。
rがそのようなgの倍数となる確率は、rがa[i]のすべてで割りきれてb[i]のいずれでも割れない確率なので1/a[i]の積に(1-(1/b[i]))の積を掛ければよい。

問題はg==1の時にcやhがどうなるかである。
フェルマーの小定理よりr^Q = 1 (mod P)なので、(r^q)を2の累乗乗するといずれ1になる。
その際、最大でもe回2乗すれば1になる。
各rに対し、何回2乗すればよいかを考えると、2^eの巡回群より以下の確率になることがわかる。

  • i=0 : 1/(2^e)
  • i=1～e ： 1/(2^(e+1+i))

hが同様にa[i]の積の倍数であってb[i]の倍数でない確率は以下のように求められる。
a[i]-1=2^e[i]*q[i]、b[i]-1=2^g[i]*s[i]と表現できるとする。(q,sは奇数)
各kにおいて以下の積：

  • 各a[i]においてk回累乗すると1になる確率の積
  • 各b[i]においてk回未満累乗した時点で1になる確率の積

g==1の時の確率に、上記積のk毎の和を掛けると、3-3の手順でhがa[i]とb[i]に分割できる確率が求められる。

int T;
int N;
ll P[20];
int bi[20];
int minbi[1<<15];
double bip[20][65];
double bipsum[20][65];
double memo[1<<15];

double dodo(int pm) {
	
	if(__builtin_popcount(pm)<=1) return __builtin_popcount(pm);
	if(memo[pm]>=0) return memo[pm];
	
	vector<int> id;
	double ret=0;
	double left=1, p2b=1;
	int mask,i,x;
	
	FOR(i,N) if(pm&(1<<i)) id.push_back(i), p2b *= (P[i]-1.0)/P[i];
	int n=id.size();
	
	
	FOR(mask,1<<n) {
		double p1=1;
		int nm[2]={};
		FOR(i,n) {
			if(mask&(1<<i)) p1 *= 1.0/P[id[i]], nm[0] |= 1<<id[i];
			else p1 *= (P[id[i]]-1.0)/P[id[i]], nm[1] |= 1<<id[i];
		}
		if(nm[0] && nm[1]) {
			int mb=minbi[nm[0]];
			
			double p2=0;
			for(x=1;x<=63;x++) {
				double hoge=1;
				FOR(i,n) {
					if(mask&(1<<i)) hoge *= bip[id[i]][x];
					else hoge *= bipsum[id[i]][x-1];
				}
				p2+=hoge;
			}
			
			double pat = p1 + p2*p2b;
			left -= pat;
			ret += pat * (dodo(nm[0])+dodo(nm[1]));
		}
	}
	
	return memo[pm]=(1+ret)/(1-left);
}

void solve() {
	int i,j,k,l,r,x,y; string s;
	
	cin>>N;
	FOR(i,1<<N) memo[i]=-1;
	
	ZERO(bip);
	ZERO(bipsum);
	FOR(i,N) {
		cin>>P[i];
		ll p=P[i]-1;
		bi[i]=0;
		while(p%2==0) bi[i]++, p/=2;
		
		bipsum[i][0]=bip[i][0]=pow(0.5,bi[i]);
		for(x=1;x<=64;x++) {
			if(x<=bi[i]) bip[i][x]=pow(0.5,bi[i]+1-x);
			bipsum[i][x]=bipsum[i][x-1]+bip[i][x];
		}
		
	}
	
	FOR(i,1<<N) {
		minbi[i]=100;
		FOR(j,N) if(i & (1<<j)) minbi[i]=min(minbi[i],bi[j]);
	}
	
	_P("%.12lf\n",dodo((1<<N)-1));
}

まとめ

Euler勢怖い。