AtCoder ABC #271 (京セラプログラミングコンテスト2022) : G - Access Counter

だいぶてこずった。
https://atcoder.jp/contests/abc271/tasks/abc271_g

問題

1日24時間のうち、毎時0分にウェブサイトに以下のアクセスが発生する。

• 事前に、毎時高橋君と青木君どちらがアクセスし得るか確定している。
• 高橋君がアクセスする可能性がある時間では、確率Xで実際にアクセスする
• 青木君がアクセスする可能性がある時間では、確率Yで実際にアクセスする

今0時直前時点でアクセス数が0の場合、そこからN回目のアクセスが青木君となる確率を求めよ。

解法

f(n,h1,h2) := 今h1時直前であるとき、あとn回後のアクセスがh2時のカウンタである確率
とする。
その時、青木君がアクセスする可能性がある時間帯a1,a2,a3....において、解はf(n,0,a1)+f(n,0,a2)+....となる。

• h2時に高橋君がアクセスする可能性がある場合
  • f(n,h1,h2+1) = X*f(n,h1+1,h2+1)+(1-X)*f(n-1,h1,h2+1)
• h2時に青木君君がアクセスする可能性がある場合
  • f(n,h1,h2+1) = Y*f(n,h1+1,h2+1)+(1-Y)*f(n-1,h1,h2+1)

f(n,*,*)からf(n-1,*,*)を計算する式を、行列Aの乗算で表現すると、(A^-1)^Nを求めればf(N,0,*)を計算できる。

ll N,X,Y;
string C;

const int MAT=48;
ll mo=998244353;
ll ma[MAT][MAT],ma2[MAT][MAT],V[MAT],R[MAT];
struct Mat { ll v[MAT][MAT]; Mat(){ZERO(v);};};

ll rev(ll a, ll n = mo-2) {
	ll r=1;
	while(n) r=r*((n%2)?a:1)%mo,a=a*a%mo,n>>=1;
	return r;
}

ll modpow(ll a, ll n = mo-2) {
	ll r=1;a%=mo;
	while(n) r=r*((n%2)?a:1)%mo,a=a*a%mo,n>>=1;
	return r;
}

Mat mulmat(Mat& a,Mat& b,int n=MAT) {
	ll mo2=4*mo*mo;
	int x,y,z; Mat r;
	FOR(x,n) FOR(y,n) r.v[x][y]=0;
	FOR(x,n) FOR(z,n) FOR(y,n) {
		r.v[x][y] += a.v[x][z]*b.v[z][y];
		if(r.v[x][y]>mo2) r.v[x][y] -= mo2;
	}
	FOR(x,n) FOR(y,n) r.v[x][y]%=mo;
	return r;
}

Mat powmat(ll p,Mat a,int n=MAT) {
	int i,x,y; Mat r;
	FOR(x,n) FOR(y,n) r.v[x][y]=0;
	FOR(i,n) r.v[i][i]=1;
	while(p) {
		if(p%2) r=mulmat(r,a,n);
		a=mulmat(a,a,n);
		p>>=1;
	}
	return r;
}

//行と列を決めての
int mat_rank(int r,int c,ll ma_[MAT][MAT]) {
	static ll tmp[MAT][MAT];
	int y,x,step;
	FOR(y,r) FOR(x,c) tmp[y][x]=ma_[y][x];
	int ret=0;
	FOR(step,c) {
		for(y=ret;y<r;y++) if(tmp[y][step]) break;
		if(y==r) continue;
		for(x=step;x<c;x++) swap(tmp[ret][x],tmp[y][x]);
		int re=rev(tmp[ret][step]);
		for(x=step;x<c;x++) tmp[ret][x]=tmp[ret][x]*re%mo;
		
		FOR(y,r) if(y!=ret && tmp[y][step]) {
			re=tmp[y][step];
			for(x=step;x<c;x++) tmp[y][x]=((tmp[y][x]-re*tmp[ret][x])%mo+mo)%mo;
		}
		ret++;
	}
	FOR(y,r) FOR(x,c) ma_[y][x]=tmp[y][x];
	return ret;
	
}

void solve() {
	int i,j,k,l,r,x,y; string s;
	
	cin>>N>>X>>Y>>C;
	X=X*modpow(100)%mo;
	Y=Y*modpow(100)%mo;
	FORR(c,C) c=c=='A';
	FOR(i,24) {
		ma[i][i+24]=1;
		if(C[i]==0) {
			ma[(i+1)%24][i]=modpow(X);
			ma[(i+1)%24][(i+1)%24]=(X-1)*modpow(X)%mo;
		}
		else {
			ma[(i+1)%24][i]=modpow(Y);
			ma[(i+1)%24][(i+1)%24]=(Y-1)*modpow(Y)%mo;
		}
	}
	mat_rank(24,48,ma);
	Mat ma2;
	FOR(y,24) {
		FOR(x,24) ma2.v[y][x]=ma[y][x+24];
	}
	ma2=powmat(N,ma2);
	ll ret=0;
	FOR(i,24) if(C[(i+23)%24]==1) ret+=ma2.v[0][i];
	cout<<ret%mo<<endl;
	
}

まとめ

逆行列求めるコード久々に使った。