SRM613 Med

問題:RandomGCD
概要 low以上high以下の整数をn個並べてできる(high-low+1)^n個の数列のうち,gcdがkであるものの個数をmod1e9+7で求めよ
1<=n,k,low,high<=1e9,low<=high,high-low<=1e5

自分の解法:
さすがに包除原理やるだけ…。まず全体をkで割るみたいなことをしておく.後は包除原理をするだけだが、全てについて真面目にやってると間に合わないので、10^5まではメビウス関数を計算しておき(もちろん素数でループをまわす!)、10^5より大きいやつdに関しては、高々数列は1通り(全部d)となる.
よって今度はlowからhighまでの数xに対し、まだ考慮されていない(つまり10^5より大きい)xに関連する部分のメビウス関数を計算する。(具体的にはxの素因数を列挙(高々10個)して、(2^個数)個の約数だけ確かめれば良い

公式解法:f(low,high)="low~highでgcdが1になる組み合わせ"として、メモ化再帰.

ミス:
10^5より大きいやつになりうることに気づかない
素因数列挙時に10^5より大きい素因数を忘れる

メモ:
a以上の最小のbの倍数は(a+b-1)/b*b
a以下の最大のbの倍数はa/b*b

my solution

#include <iostream>
#include <vector>
#define rep(i,n) for(int i=0;i<n;i++)
#define rep1(i,n) for(int i=1;i<=n;i++)
#define all(c) c.begin(),c.end()
#define pb push_back
#define fs first
#define sc second
using namespace std;
typedef long long ll;
ll mod=1e9+7;
const ll ccc=100000;
bool prime[ccc+1];
vector<ll> pr;
int t[ccc+1];
vector<int> d,facd[ccc+1];
void makeprime(){
	ll i,j;
	for(i=2;i<=ccc;i++) prime[i]=true;
	for(i=2;i*i<=ccc;i++) if(prime[i]) for(j=2;j*i<=ccc;j++) prime[j*i]=false;
	for(i=2;i<=ccc;i++) if(prime[i]) pr.push_back(i);
}
int tmp[ccc+1];
class RandomGCD{
	public:
	int countTuples(int n, int k, int low, int high){
		low=(low+k-1)/k,high=high/k;
		if(high-low<0) return 0;
		makeprime();
		rep(i,ccc+1) t[i]=1;
		rep(i,pr.size()){
			int p=pr[i];
			for(int j=1;j*p<=ccc;j++){
				if(j%p==0) t[j*p]=0;
				else t[j*p]*=-1;
			}
		}
		ll ans=0;
		for(int i=1;i<=ccc;i++){
			ll x=high/i-(low-1)/i,now=1;
			for(int j=0;(n>>j)>0;j++){
				if((n>>j)&1) now=now*x%mod;
				x=x*x%mod;
			}
			ans+=t[i]*now;
		}
		rep(i,high-low+1) tmp[i]=i+low;
		rep(i,pr.size()){
			int p=pr[i];
			for(int j=(low+p-1)/p;j*p<=high;j++){
				facd[j*p-low].push_back(p);
				while(tmp[j*p-low]%p==0) tmp[j*p-low]/=p;
			}
		}
		rep(i,high-low+1) if(tmp[i]>1) facd[i].push_back(tmp[i]);
		rep(i,high-low+1){
			int s=facd[i].size();
			rep(j,1<<s){
				ll now=1;
				int cnt=1;
				rep(h,s) if((j>>h)&1) now*=facd[i][h],cnt=-cnt;
				if(now>ccc) ans+=cnt;
			}
		}
		return (ans%mod+mod)%mod;
	}
};