素因数分解 - yanprog @ ウィキ

単独の整数(n<2^31)を素因数分解する。使う場面は多分あまりない。それなりに早い。
factor.is_prime(ll n) ｎが素数かどうか判定。
prime_divisor(ll n) nを素因数分解(ソートして返す)。例えばn=60なら[2,2,3,5]
factor.phi(ll n) ｎのオイラー関数
factor.divisor(ll n) nの約数全体を返す。例えば n=2なら [1,2,5,10]

struct FACTOR{
	static const int table_size=65536;
	vvl table;
	FACTOR(){
		table.resize(table_size);
		for(int p=2; p<table_size; p++){
			if(table[p].size()>0)continue;
			for(int k=p; k<table_size; k+=p){
				int k2=k;
				while(k2%p==0){
					table[k].pb(p);
					k2/=p;
				}
			}
		}
	}
 
	ll pow(ll a, ll n, ll m){
		if(n==0)return 1;
		return pow(a*a%m, n/2, m)  * (n%2==0 ? 1 : a ) % m;
	}
 
	unsigned int gcd(unsigned int a, unsigned int b){
		return a==0 ? b : gcd(b%a, a);
	}
 
	bool miller_rabin_prime(ll n, ll a){
		if(n%a==0)return false;
		if(pow(a,n-1,n)!=1)return false;
		ll s=1, d=n-1;
		while(d%2==0){
			d/=2;
			s+=1;
		}
		ll v=pow(a,d,n);
		if(v==1)return true;
		while(true){
			if(v+1==n)return true;
			if(v==1)return false;
			v=v*v%n;
		}
	}
 
	bool is_prime(ll n){
		assert(n>0);
		if(n<table_size)return table[n].size()==1;
		if(miller_rabin_prime(n,2) && miller_rabin_prime(n,7) && miller_rabin_prime(n,61))return true;
	}
 
	void prime_divisor(ll n, vl &ret){
		while(n%2==0){
			ret.pb(2);
			n/=2;
		}
		while(n%3==0){
			ret.pb(3);
			n/=3;
		}
		while(n%5==0){
			ret.pb(5);
			n/=5;
		}
		prime_divisor_after(n,ret);
	}
 
	void prime_divisor_after(ll n, vl &ret){
		assert(n>0);
		if(n<table_size){
			REP(i,table[n].size())ret.pb(table[n][i]);
			return;
		}
		if(is_prime(n)){
			ret.pb(n);
			return;
		}
		ll a=1,b=a*a+1;
		int step=0;
		while(true){
			ll d=gcd(abs(a-b), n);
			if(d!=1 && d!=n){
				prime_divisor_after(d, ret);
				prime_divisor_after(n/d, ret);
				return;
			}
 
			b=(b*b+1)%n;
			b=(b*b+1)%n;
			a=(a*a+1)%n;
 
			step++;
			if(a==b){
				a=step%n;
				b=(a*a+1)%n;
			}
		}
	}
 
	vl divisor(ll n){
		vl ps;
		prime_divisor(n, ps);
		map<ll, int> cnt;
		REP(i,ps.size()) cnt[ps[i]]++;
		vl ret;
		ret.pb(1);
		EACH(pp, cnt){
			vl ds;
			ll v=1;
			for(int k=0; k<=pp->second; k++){
				ds.pb(v);
				v*=pp->first;
			}
			vl ret2;
			REP(i1, ret.size()) REP(i2, ds.size()) ret2.pb(ret[i1]*ds[i2]);
			ret=ret2;
		}
		sort(ret.begin(), ret.end());
		return ret;
	}
 
	ll phi(ll n){
		vl ps;
		prime_divisor(n, ps);
		set<ll> set_ps(ps.begin(), ps.end());
		EACH(pp, set_ps){
			ll p=*pp;
			n=n/p*(p-1);
		}
		return n;
	}
 
} factor;
vl prime_divisor(ll n){
	vl ret;
	factor.prime_divisor(n,ret);
	sort(ret.begin(), ret.end());
	return ret;
}