これ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