結果

問題 No.75 回数の期待値の問題
コンテスト
ユーザー koyumeishi
提出日時 2015-03-05 00:35:13
言語 C++11（廃止可能性あり）
(gcc 13.3.0)
結果
AC  
実行時間 13 ms / 5,000 ms
コード長 3,473 bytes
コンパイル時間 957 ms
コンパイル使用メモリ 84,308 KB
実行使用メモリ 6,944 KB
最終ジャッジ日時 2024-06-24 07:59:35
合計ジャッジ時間 1,647 ms
ジャッジサーバーID
（参考情報） 		judge3 / judge5
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
sample AC * 4
other AC * 16
権限があれば一括ダウンロードができます

ソースコード

diff #
raw source code 

#include <iostream>
#include <vector>
#include <cstdio>
#include <sstream>
#include <map>
#include <string>
#include <algorithm>
#include <queue>
#include <cmath>
#include <set>
using namespace std;

//[n*p] * [p*m] => [n*m]
template<class T>
vector< vector<T> > multmat(const vector<vector<T> > &A, const vector<vector<T>> &B, int n, int p, int m){
	vector<vector<T> > C(n, vector<T>(m,0));
	for(int i=0; i<n; i++){
		for(int k=0; k<p; k++){
			for(int j=0; j<m; j++){
				C[i][j] += A[i][k] * B[k][j];
				//C[i][j] %= mod;
			}
		}
	}
	return C;
}

//A[n*n]^k 
template<class T>
vector< vector<T> > mat_pow(vector<vector<T> > A, int k){
	int n = A.size();
	vector<vector<T> > ret(n, vector<T>(n, 0) );
	for(int i=0; i<n; i++){
		ret[i][i] = 1;
	}
	while(k>0){
		if(k&1) ret = multmat(A,ret, n,n,n);
		A = multmat(A,A, n,n,n);
		k>>=1;
	}
	return ret;
}


//A[n*p]*b[p*m] = x[n*m]
vector<vector<double> > gaussian_elimination(const vector<vector<double> > &A, vector<vector<double> > &x){
	bool valid = true;
	int rank=0;
	int n = A.size();
	int p = A[0].size();
	vector< vector<double> > R = A;
	for(int i=0; i<n; i++){
		R[i].insert(R[i].end(), x[i].begin(), x[i].end());
	}
	int m = R[0].size();	//p + x[0].size()

	//foward
	//row
	for(int i=0; i<n; i++){
		//pivot
		int pivot_row = i;
		double pivot_val = fabs(R[i][i]);
		//choose pivot row
		for(int j=i+1; j<n; j++){
			if(pivot_val < fabs(R[j][i])){
				pivot_row = j;
				pivot_val = fabs(R[j][i]);
			}
		}
		if(pivot_row != i){
			swap(R[i], R[pivot_row]);
		}

		if(pivot_val <= 1e-9){
			
			continue;
		}

		rank++;

		{
			double divisor = 1.0/R[i][i];
			//A[i][j] == 0, such that (0 <= j < i)
			for(int j=i; j<m; j++){
				R[i][j] *= divisor;
			}
		}
		
		for(int j=i+1; j<n; j++){
			double divisor = R[j][i];	// R[j][i]/R[i][i], but R[i][i] is 1.0
			for(int k=i; k<m; k++){
				R[j][k] -= R[i][k] * divisor;
			}
		}
	}

	//backward
	//row
	for(int i=n-1; i>=0; i--){
		for(int j=i-1; j>=0; j--){
			double c = R[j][i];
			for(int k=i; k<m; k++){
				R[j][k] -= R[i][k] * c;
			}
		}
	}

	for(int i=0; i<n; i++){
		bool row_valid = false;
		for(int j=i; j<m; j++){
			if(fabs(R[i][j]) <= 1e-9){
				R[i][j] = 0.0;
			}else if(j<p){
				row_valid = true;
			}
		}
		if(row_valid == false && R[i][n] != 0.0){
			valid = false;
		}
	}

	return R;
}


//#define DBG_

void print_vec(const vector<vector<double>> &v){
#ifdef DBG_
	for(auto l : v){
		for(auto x : l){
			cerr << x << "\t";
		}
		cerr << endl;
	}
	cerr << endl;
#endif
}

int main(){
	int K;
	cin >> K;

	double p = 1.0/6.0;
	vector<vector<double>> P(K+1, vector<double>(K+1, 0));
	for(int i=0; i<K; i++){
		for(int k=1; k<=6; k++){
			if(i+k<=K) P[i][i+k] = p;
			else P[i][0] += p;
		}	
	}
	P[K][K] = 1;

	//cerr << "P" << endl;
	print_vec(P);

	vector<vector<double>> E(K, vector<double>(K, 0));
	for(int i=0; i<K; i++) E[i][i] = 1;

	//cerr << "E" << endl;
	print_vec(E);

	auto Q = E;
	for(int i=0; i<K; i++){
		for(int j=0; j<K; j++){
			Q[i][j] -= P[i][j];
		}
	}
	//cerr << "Q" << endl;
	print_vec(Q);

	auto M_ = gaussian_elimination(Q, E);
	vector<vector<double>> M(K, vector<double>(K, 0));
	for(int i=0; i<K; i++){
		for(int j=0; j<K; j++){
			M[i][j] = M_[i][K+j];
		}
	}


	//cerr << "M" << endl;
	print_vec(M);

	vector<vector<double>> g(K, vector<double>(1,1));

	//cerr << "g" << endl;
	print_vec(g);

	auto X = multmat<double>(M,g, K,K,1);

	//cerr << "X" << endl;
	print_vec(X);

	cout << X[0][0] << endl;
	return 0;
}
0