yukicoder No.301 サイコロで確率問題 (1)

問題

No.301 サイコロで確率問題 (1) - yukicoder

解法

N<5000ならDPで計算、N≧5000ならN+5/3を出力すればいい。
Nが大きいときの振る舞いは実験すれば分かる。

考察

一次式でどうして近似できるのか真面目に考えてみた。

漸化式
{E_{n+6}=(E_{n+5}+\cdots+E_n)/6+1}
の一般解を求めることを考える。定数項が邪魔なので階差数列{E'_n}で漸化式を立てなおすと、
{E'_{n+6}=(E'_{n+5}+\cdots+E'_n)/6}
となり、定数係数の線形漸化式になる。

特性方程式を某サイトに投げて解かせてみると、
{x_1=1}
{x_2=-0.670332}
{x_3=-0.375695-0.570175i}
{x_4=-0.375695+0.570175i}
{x_5=0.294195-0.668367i}
{x_6=0.294195+0.668367i}
だと分かるので、一般解は
{E'_n=c_1x_1^n+c_2x_2^2+...+c_6x_6^n}

{x_1=1}であることに気をつけて{E_n}を復元すると
{E_n=d_1n+d_2x_2^n+...+d_6x_6^n+e}
という形になる。

絶対値を計算してみると|x2|,...,|x6|<0.74が分かるので、nが十分大きいとき{x_i^n}の項は0に限りなく近づく。

（説明を追加しました　2015-11-15 23:10）
数列{X_n},{Y_n}を次のように定義する。

{X_{i+6}=(X_{i+5}+\cdots+X_i)/6}
{Y_{i+6}=(Y_{i+5}+\cdots+Y_i)/6+1}
{X_0=0,X_i=1\;(i<0),Y_0=0,Y_i=0\;(i<0)}

このように定義すると求める期待値は
{dp[i]=X_i dp[N]+Y_i}
という形に表せる。

ここで{X_n},{Y_n}は{N}に依らない数列なので、ある定数a,b,cが存在して、{X_i}は{a}に、{Y_i}は{bi+c}に近づく。

{dp[N]=a\cdot dp[N]+bN+c}
{dp[N]=(bN+c)/(1-a)}
以上より{dp[N]}は一次式で近似できる。

ソースコード

本番中に書いた行列累乗をそのまま使って書いた。

#include <bits/stdc++.h>
#define GET_MACRO(a, b, c, NAME, ...) NAME
#define rep(...) GET_MACRO(__VA_ARGS__, rep3, rep2)(__VA_ARGS__)
#define rep2(i, a) rep3 (i, 0, a)
#define rep3(i, a, b) for (int i = (a); i < (b); i++)
#define repr(...) GET_MACRO(__VA_ARGS__, repr3, repr2)(__VA_ARGS__)
#define repr2(i, a) repr3 (i, 0, a)
#define repr3(i, a, b) for (int i = (b) - 1; i >= (a); i--)
#define chmin(a, b) ((b) < a && (a = (b), true))
#define chmax(a, b) (a < (b) && (a = (b), true))
using namespace std;
typedef long long ll;

typedef long double D;

const int maxdim = 8;
int dim;
typedef D mat[maxdim][maxdim];
mat M0, M1, M2, M3, M4;
mat Mp[80];

void matmul(mat A, mat B, mat res) {
    rep (i, dim) rep (j, dim) M0[i][j] = 0;
    rep (i, dim) rep (k, dim) rep (j, dim) {
        M0[i][j] += A[i][k] * B[k][j];
    }
    rep (i, dim) rep (j, dim) res[i][j] = M0[i][j];
}

void matpow(mat A, ll b, mat res) {
    rep (i, dim) rep (j, dim) M1[i][j] = A[i][j];
    rep (i, dim) rep (j, dim) M2[i][j] = i == j;
    int i = 1;
    while (b > 0) {
        if (b & 1) matmul(Mp[i], M2, M2);
        b >>= 1;
        i++;
    }
    rep (i, dim) rep (j, dim) res[i][j] = M2[i][j];
}
void companion(vector<D> a, mat res) {
    int n = a.size();
    rep (i, n - 1) res[i][i + 1] = 1;
    rep (i, n) res[n - 1][i] = a[i];
    res[n][1] = res[n][n] = 1;
}

int main() {
    int T;
    cin >> T;
    dim = 7;
    D p = (long double)1.0 / 6;
    companion(vector<D>(6, p), M4);
    rep (i, 7) rep (j, 7) Mp[1][i][j] = M4[i][j];
    rep (i, 1, 70) {
        matmul(Mp[i], Mp[i], Mp[i + 1]);
    }
    rep (i_, T) {
        ll N;
        cin >> N;
        if (N < 1000) {
            matpow(M4, N + 4, M3);
            D s = M3[6][4], t = M3[6][5];
            D ans = t / (s - t * 5.0 / 6.0);
            cout << fixed << setprecision(50) << ans << endl;
        } else {
      　　　cout << (N + 1) << ".66666666666666666666666" << endl;
        }
    }
    return 0;
}

コメント

本番中は行列累乗でE[n]を求めてたけどnが大きい時の精度がダメダメだった。
long double万能だと思ってたけど、整数部分もでかいから精度足りなくなるみたい。

実験すれば確かに気づくし実験は重要。