SRM683

Programming

わりと面白かった

反省点:
easyでオーバーフローさせた
medでauto&とすべきところをautoにしていた
hardで配列外アクセスを起こした(これはコンテスト後)



Easy:
N個の正整数が円状に並んでいる.隣り合ったものに+1,-1を足すことができる.最小で何回操作すれば望みの数列にできるか N<1000,値≤1e9


最適解(のうち一つ)では、ある辺(隣り合った二点)が存在し,その辺に関して操作をしない.(∵辺にどっち何個に移動させるかを書いとく.一周ぐるっと+1みたいなことをしてもできる数列は変わらない.もし0がなかったとしたら,辺の向いてる向きの少ない方に向かって1ずつ流していくと操作回数は減る.どこかの辺が0になるまでこれで改善できる)

したがって、その辺を全探索して,そうするとそこで円が切れるので一意に定まる O(N^2)


Med:
正整数がN個書かれている.2つの書かれている整数を選び,その2つを消して,代わりにgcdとlcmを書く という操作を好きなだけできる. 総和が最大になるようにせよ (その時の総和 mod 1e9+7を答えよ). N<1e5,書かれている数<1e7

a=2^k1 * 3^k2 * ... ,b=2^l1 * ...と素因数分解(指数0も明示的に書く)すると,lcm=2^max(k1,l1) * 3^max(k2,l2) * ..., gcd=2^min(k1,l1) * 3^min(k2,l2) * ... となる.上の式から,各素数に対して,指数のmultisetは操作のあとで変わらないことがわかる.さらに操作によってでかいものを一つに集めることが可能.lcm*gcd=a*bであるから,和を大きくするにはできるだけ大きいやつと小さい奴の差を激しくすると良い.したがって各素因数に対して指数たちをソートし,でかい順に並べてかけたものを足せば良い.1まで全部愚直にかけてるとTLEなので,ちゃんと指数が1以上のもののみを見れるようにmapとかで見ていき終わったら消すみたいな操作が必要.

Hard:
二次元平面上のある格子点(x0,y0)からKターン4方向にランダムウォークする.(x^m * y^n の期待値)*4^K を mod 1e9+7で求めよ.
K<1e9, x0,y0の絶対値<1e9,n,m<100


解法が3つぐらいあるので紹介.
まずxとyは独立でないので期待値の積は積の期待値にはならない.SRM573d1Hardの考え方のように,45度回転してX=x+yとY=x-yは独立になる.(ここまで共通)

1.
x^m * y^n をX,Yの多項式で表す. あとはX^a,Y^bの期待値がわかればすべてが求まる.X_i=(x0+a0+a1+...+a_{i-1})と書く(iターン目にX方向にいくら動いたか がa_{i}).K=iの時のX^0,X^1,..X^aの期待値(このベクトルをvx_{i}とおく)がわかっているとして,i+1の時にどうなるかというと,X_iをa_{i}の次数別に展開すると,他の部分はvx_{i-1}から計算できる.従って,vx_{i}からvx_{i-1}を計算する行列がわかる.この行列はiによらないので,行列累乗するとO(N^3 log K)で答えが求まる.適当にやるとTLEっぽいが行列が下三角であることなどを使って計算をちょっと略すと通ったりするらしい.(なんらかのsum)^N を求める時の一般的なテク.

2.
行列を求めるとこまでは1と一緒. (K回動いた時のx^m * y^nの期待値)はなんか高々m+n(?)次式らしいので,行列累乗しなくても,K=0,...m+n+1くらいでの答えを求めてO(N^3)、あとは多項式補間(O(N^2))すればOK←魔剤！？ (多項式になる証明はわかってません.(snukeは実験してみたと言っていた こういう系ではわりと一般的なテクかもしれない)).

3.
x^m * y^n を直接x0,y0,a_{i},b_{i}の式で書き下す.
x0*=2,y0*=2
(x0+a0+...+a_{K-1}+b0+..b_{K-1})^n * (y0+a0+...+a_{K-1}-b0-..b_{K-1})^m /(2^(n+m))
これのa_{i},b_{i}にそれぞれ+-1を入れた4^K通りの値の総和を求めればよい.ここで式を展開した時にa_{i}やb_{i}の奇数乗の項があると,+1と-1を入れた時に打ち消し合うため,全部偶数乗の項しか考えなくていい.このような項では値はa,bに関する部分の値は全部1なので,そのような項の個数だけ求めれば良い,といいたいが,右の方に-b_{i}みたいなのがあるのでそこの偶奇には気をつけないといけない.
どうやるかというと,まずx0^(N-i) とy0^(M-j)のi,jを固定する(これしないと流石に総和が出せそうにない(N,Mから引いてるのは後の都合)).そして,^nの方でa_{}から奇数次のものがk個,b_{}からl個 のk,lを固定する.この時^mの方でもおなじになってることがすべて偶数次の必要十分条件.これで全体の符号がlの偶奇で判別できる.この項に関していくら足せばいいかというと,
N_C_(N-i) * M_C_(M-j) * x0^i * y0^j * calc(i,k,l) * calc(j,k,l) / K_C_k / K_C_l
です.
左4つは二項定理のとこで,後半は,まず
calc(n,i,j)=(a_0+..a_{K-1}+b_0+..b_{K-1})^n で「aの方に奇数次のものがi個,bの方に奇数次のものがj個ある」ような項の数.
K_C_kとかで割ってる理由は,calcはどのi個,どのj個が奇数次か,を無視していて,それが異なったら全体が偶数次にならないため,一致する分だけを計算するようにする係数.

calcはどう計算すればいいかというと,普通にDP O(N^3) n,i,jからの遷移は,(a,b)の方から今(奇数,偶数)次のものを取る を全部試せばO(1).

全体ではi,j,k,lのループを回してるのでO(N^4)で解けた.

なるほどなあ(3はrngさんの解法)


問題を解くにはいろいろな方法があるなあ(でもやっぱり基本はわからなかったらちゃんと数式にすることだと思う)

3で解きました

#include <bits/stdc++.h>
#define rep(i,n) for(int i=0;i<(int)(n);i++)
#define rep1(i,n) for(int i=1;i<=(int)(n);i++)
#define all(c) c.begin(),c.end()
#define pb push_back
#define fs first
#define sc second
#define show(x) cout << #x << " = " << x << endl
#define chmin(x,y) x=min(x,y)
#define chmax(x,y) x=max(x,y)
using namespace std;
typedef long long ll;
ll mod=1e9+7;
ll px[101],py[101];
ll way[101][101][101];
ll C[101][101];
ll iCT[101];
void add(ll &x,ll y){
	x+=y;
	x%=mod;
}
ll ex(ll x,ll p){
	ll a=1;
	while(p){
		if(p%2) a=a*x%mod;
		x=x*x%mod;
		p/=2;
	}
	return a;
}
ll inv(ll x){
	return ex(x,mod-2);
}
class RandomWalkOnGrid{
	public:
	int getExpectation(int x0, int y0, int T, int N, int M){
		x0=x0*2%mod,y0=y0*2%mod;
		px[0]=1,py[0]=1;
		rep(i,100) px[i+1]=px[i]*x0%mod,py[i+1]=py[i]*y0%mod;
		rep(i,101){
			C[i][0]=C[i][i]=1;
			for(int j=1;j<i;j++) C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod;
		}
		iCT[0]=1;
		rep(i,min(100,T)){
			iCT[i+1]=iCT[i]*(i+1)%mod*inv(T-i)%mod;
		}
		way[0][0][0]=1;
		rep(n,max(N,M)){
			rep(a,n+1) rep(b,n+1-a){
				if(a>T||b>T) continue;
				ll v=way[n][a][b];
				if(v==0) continue;
				if(a<T) add(way[n+1][a+1][b],v*(T-a));
				if(b<T) add(way[n+1][a][b+1],v*(T-b));
				if(a>0) add(way[n+1][a-1][b],v*a);
				if(b>0) add(way[n+1][a][b-1],v*b);
			}
		}
		ll ans=0;
		rep(i,N+1) rep(j,M+1){		//x0^(N-i) * y0^(M-j)
			int mn=min(i,j);
			rep(k,mn+1) rep(l,mn+1-k){
				if(k>T||l>T) continue;
				ll tmp=px[N-i]*py[M-j]%mod*C[N][i]%mod*C[M][j]%mod*way[i][k][l]%mod*way[j][k][l]%mod*iCT[k]%mod*iCT[l]%mod;
				if(l%2) tmp=-tmp;
				add(ans,tmp);
			}
		}
		ll i2=(mod+1)/2;
		rep(i,N+M) ans=ans*i2%mod;
		ans=ans*ex(4,T);
		return (ans%mod+mod)%mod;
	}
};