概率dp Lights inside 3D Grid

Lights inside 3D Grid
You are given a 3D grid, which has dimensions X, Y and Z. Each of the X x Y x Z cells contains a light. Initially all lights are off. You will have K turns. In each of the K turns,
1.You select a cell A randomly from the grid,
2.You select a cell B randomly from the grid and
3.Toggle the states of all the bulbs bounded by cell A and cell B, i.e. make all the ON lights OFF and make all the OFF lights ON which are bounded by A and B. To be clear, consider cell A is (x1, y1, z1) and cell B is (x2, y2, z2). Then you have to toggle all the bulbs in grid cell (x, y, z) where min(x1, x2) ≤ x ≤ max(x1, x2), min(y1, y2) ≤ y ≤ max(y1, y2) and min(z1, z2) ≤ z ≤ max(z1, z2).

Your task is to find the expected number of lights to be ON after K turns.

Input

Input starts with an integer T (≤ 50), denoting the number of test cases.
Each case starts with a line containing four integers X, Y, Z (1 ≤ X, Y, Z ≤ 100) and K (0 ≤ K ≤ 10000).

Output

For each case, print the case number and the expected number of lights that are ON after K turns. Errors less than 10-6 will be ignored.

Sample Input

5
1 2 3 5
1 1 1 1
1 2 3 0
2 3 4 1
2 3 4 2

Sample Output

Case 1: 2.9998713992
Case 2: 1
Case 3: 0
Case 4: 6.375
Case 5: 9.09765625

遍历每个点v(x,y,z)，对P的每维，若x轴是a这个数，

a被取到的概率$p_x=1-\frac{(X-a{)}^2+(a-1{)}^2}{X^2}$

执行一次打开或关闭点v的概率$p=p_x{p}_y{p}_z$

$v_{ans}=p^1(1-p{)}^k-1+p^3(p-1{)}^{k-3}+\cdots +p^m(1-p{)}^{k-m}$，其中，m是

小于k的最大奇数。

设$f(k)=v_{ans}$：k轮里，点亮次数为奇数次的概率

$g(k)$:k轮里，点亮次数为偶数次的概率，所以$g(k)+f(k)=1$

又有$f(k)=g(k-1)p+f(k-1)(1-p)$

所以$f(k)=(1-2p)f(k-1)+p$
能推出通项公式，也能直接矩阵。

#include<bits/stdc++.h>
using namespace std;
#define Init(a,v) memset(a,v,sizeof(a))
#define Tie ios::sync_with_stdio(0);cin.tie(0)
typedef long long ll;
const ll MAXN=1e1+8,inf=0x3f3f3f3f,mod=1e9+7;
int a[MAXN];
const int SIZE=2;
struct matrix{
    double a[SIZE][SIZE];
    matrix(){memset(a,0,sizeof(a));}
    matrix operator*(matrix&o){
        matrix res;
        for(int i=0;i<SIZE;++i)for(int j=0;j<SIZE;++j)
        for(int k=0;k<SIZE;++k)res.a[i][j]+=a[i][k]*o.a[k][j];
        return res;
    }
    matrix pow(ll k){
        matrix res,x=*(this);
        for(int i=0;i<SIZE;++i)res.a[i][i]=1;
        while(k){
            if(k&1)res=res*x;
            x=x*x;k>>=1;
        }return res;
    }
};
int C[105];
inline double f(int x,int a){
    double p=1.0*(C[a-1]+C[x-a])/C[x];
    return 1-p;
}
int main() {
    Tie;
    for(int i=1;i<=100;++i)C[i]=i*i;
    int t,cas=1;
    cin>>t;
    int x,y,z,k;
    while(t--){
        cin>>x>>y>>z>>k;
        if(k==0){printf("Case %d: 0\n",cas++);continue;}
        double ans=0;
        for(int i=1;i<=x;++i){
            double px=f(x,i);
            for(int j=1;j<=y;++j){
                double py=f(y,j);
                for(int l=1;l<=z;++l){
                    double p=px*py*f(z,l);
                    matrix e;
                    e.a[0][0]=1-2*p;
                    e.a[0][1]=p;
                    e.a[1][1]=1;
                    e=e.pow(k-1);
                    ans+=p*e.a[0][0]+e.a[0][1];
                }
            }
        }
        printf("Case %d: %lf\n",cas++,ans);
    }
    return 0;
}