kmjp's blog

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

AtCoder ARC #123 : F - Insert Addition

これ1000ptより難しいと思うんだけど、AGCの1000ptとARCの1000ptで基準が違うのかな。
https://atcoder.jp/contests/arc123/tasks/arc123_f

問題

数列Pに対し、f(P)とはPの要素間に両隣の要素の和となる値を挿入した数列を示す。
初期状態の数列A=(a,b)と、整数Nが与えられる。
Aをf(A)に置き換えることをN回行ったうえで、AからNより大きな値を取り除く。
区間[L,R]が指定されるので、A[L]~A[R]を求めよ。

解法

A=f(A)と置き換えることを繰り返すと、いずれ追加される値が大きくなり、最終的に追加される値が(N+1)以上になるので、「N回行う」は余り意味がない。

先にAが最終的にどうなるかを示すと、
A=*1とAの各要素を2値のペアで表したとき、その間に両隣のペアの要素ごとの和が挿入される。
ペア(x,y)に対し、元の値としてはax+byが相当する。

この時、Aに入る値は、以下のようになる。

  • ax+by≦Nとなる(x,y)のうち、x,yが互いに素なものが列挙される。
  • y/xが昇順に並ぶ。

後者がポイントで、Aのn番目を求めたい場合、y/xの値を二分探索すればよい。
y/x=cを定めたとき、ax+by≦Nかつy/x≦cかつ(x,y)が互いに素なケースを列挙することになるが、これはメビウスの反転公式を使いO(N)で列挙できる。

とはいえ、これだとAの1要素を求めるのにO(NlogN)かかるので、要素毎にこれを行うとO((R-L)NlogN)かかり間に合わない。
そこで、y/xの値について、L番目の値とR番目の値を求め、あとは各xに対し、yがその範囲に収まるものを列挙し、y/xの順でソートした。

ll A,B,N;
ll L,R;

ll F[603030];
ll G[603030];

const int prime_max = 300000;
int MU[prime_max+1];
int num[prime_max+1];
int Yma[300005];
vector<int> MUs[prime_max+1];

void mebius() {
	int i,j;
	for(i=1;i<=prime_max;i++) MU[i]=1, num[i]=i;
	for(int i=2;i<=prime_max;i++) if(num[i]==i) {
		for(j=i;j<=prime_max;j+=i) {
			int x=0;
			MU[j]=-MU[j];
			while(num[j]%i==0) {
				x++;
				num[j]/=i;
			}
			if(x>=2) MU[j]=0;
		}
	}
	for(i=1;i<=prime_max;i++) if(MU[j]) {
		for(j=i;j<=prime_max;j+=i) {
			if(MU[j]>0) MUs[j].push_back(i);
			else MUs[j].push_back(-i);
		}
	}
}

ll count(double c) {
	ZERO(F);
	int x,y,i;
	for(x=0;Yma[x]>=0;x++) {
		int y=min((ll)Yma[x],(c>=1e8)?N:(ll)floor(c*x));
		F[A*x]++;
		F[A*x+(y+1)*B]--;
	}
	
	for(i=B;i<=N;i++) F[i]+=F[i-B];
	F[0]--;
	//FOR(i,N+1) cout<<F[i]<<" ";
	for(i=1;i<=N;i++) F[i]+=F[i-1];
	ll ret=0;
	for(i=1;i<=N;i++) {
		if(MU[i]==1) ret+=F[N/i];
		else if(MU[i]==-1) ret-=F[N/i];
	}
	
	//cout<<"!"<<c<<" "<<ret<<endl;
	return ret;
	
}

pair<ll,ll> hoge(ll V) { // V番目のfarey seq
	double L=0,R=1e9;
	int i;
	FOR(i,100) {
		double m=(L+R)/2;
		ll cnt=count(m);
		if(cnt>=V) R=m;
		else L=m;
	}
	
	ll X=0,Y=0;
	double dif=1;
	if(R>1e7) return {0,1};
	for(i=1;i*A<=N;i++) {
		double a=round(i*R);
		if(abs(a/i-R)<dif) {
			dif=abs(a/i-R);
			X=i;
			Y=a;
		}
	}
	//cout<<V<<" "<<R<<" "<<X<<" "<<Y<<endl;
	return {X,Y};
}


void solve() {
	int i,j,k,l,r,x,y; string s;
	
	mebius();
	
	cin>>A>>B>>N>>L>>R;
	
	MINUS(Yma);
	for(i=0;1LL*i*A<=N;i++) Yma[i]=(N-i*A)/B;
	
	
	pair<ll,ll> X=hoge(L);
	if(L==R) {
		cout<<A*X.first+B*X.second<<endl;
		return;
	}
	pair<ll,ll> Y=hoge(R);
	assert(X!=Y);
	
	vector<pair<double,ll>> V;
	if(Y.first==0) V.push_back({1e9,(0<<20)+1});
	double ax=X.second*1.0/X.first;
	double bx=(Y.first==0?1e9:Y.second*1.0/Y.first);
	for(i=1;i*A<=N;i++) {
		ll ma=min((ll)floor(i*bx+1e-9),(ll)Yma[i]);
		ll mi=ceil(i*ax-1e-9);
		for(x=mi;x<=ma;x++) {
			if(__gcd(x,i)==1) V.push_back({x*1.0/i,((ll)i<<20)+x});
		}
	}
	
	assert(V.size()==R-L+1);
	sort(ALL(V));
	FORR2(a,b,V) {
		x=b>>20;
		y=b%(1<<20);
		cout<<A*x+B*y<<" ";
	}
	
	
	/*
	cout<<A*X.first+B*X.second<<" ";
	cout<<A*Y.first+B*Y.second<<" ";
	for(ll a=L+2;a<=R;a++) {
		ll p=(N+X.second)/Y.second*Y.first-X.first;
		ll q=(N+X.second)/Y.second*Y.second-X.second;
		
		if(A*p+B*q<=N) {
			cout<<A*p+B*q<<" ";
		}
		else {
			a--;
		}
		
		X=Y;
		Y={p,q};
	}
	*/
	cout<<endl;
	
	
}

まとめ

各パート知識不足で難しい。1000ptには思えない…。

*1:1,0),(0,1