CSU 1848：3-sided dice 高斯消元(2010 Southwestern European Regional Contest)

1848: 3-sided dice Submit Page    Summary    Time Limit: 1 Sec     Memory Limit: 128 Mb     Submitted: 42     Solved: 5    

---

Description

　　Just like every fall, the organizers of the Southwestern Europe Dice Simulation Contest are busy again this year. In this edition you have to simulate a 3-sided die that outputs each of three possible outcomes (which will be denoted by 1; 2 and 3) with a given probability, using three dice in a given set. The simulation is performed this way: you choose one of the given dice at random, roll it, and report its outcome. You are free to choose the probabilities of rolling each of the given dice, as long as each probability is strictly greater than zero. Before distributing the mate- rials to the contestants, the organizers have to verify that it is actually possible to solve this task.

　　For example, in the first test case of the sample input you have to simulate a die that yields outcome 1; 2 and 3 with probabilities 310,410 310,410 and 310 310 .We give you three dice, and in this case the i-th of them always yields outcome i, for each i = 1,2,3. Then it is possible to simulate the given die in the following fashion: roll the first die with probability 310 310 .the second one with probability 410 410 and the last one with probability 310 310 .

Input

　　The input consists of several test cases, separated by single blank lines. Each test case consists of four lines: the first three of them describe the three dice you are given and the last one describes the die you have to simulate. Each of the four lines contains 3 space-separated integers between 0 and 10 000 inclusive. These numbers will add up to 10 000, and represent 10 000 times the probability that rolling the die described in that line yields outcome 1, 2 and 3, respectively.

　　The test cases will finish with a line containing only the number zero repeated three times (also preceded with a blank line).

Output

　　For each case, your program should output a line with the word YES if it is feasible to produce the desired　die from the given ones, and NO otherwise.

Sample Input

0 0 10000
0 10000 0
10000 0 0
3000 4000 3000

0 0 10000
0 10000 0
3000 4000 3000
10000 0 0

0 0 0

Sample Output

YES
NO

        借用这次机会学习了一下高斯消元，其实呢在现场我并不会高斯消元，只是拿者队友的模板直接敲的，所以WA了也没什么奇怪的。

        高斯消元呢其实就是线性代数里的矩阵的各种初等变换，通过这些把一个增广矩阵转化为上三角矩阵，然后逐个方程的求解迭代以求出最终的解。当然，在接的过程中，由于矩阵可能非满秩，所以可能会出现有自由变元，高斯消元在迭代求解的时候会判断当前方程未知数的个数，若个数大于1，则说明存在自由变元，并会标记自由变元。当然，如果出现某个方程左边全为0但是右边不为0，自然是无解。我们规定，无解时返回-1，有唯一整数解的时候返回0，唯一浮点数解的返回-2，有无穷解的返回自由变元的个数，即基础解系的个数，也即n-rank(A)。

        然后，我们针对这题说说吧。当时跟队的学长推荐我们先做这题，然后就(￣▽￣)"呵呵了……

        题目非常简单，有三个只有三个面的骰子(别管为什么只有三个面)，然后告诉你每个分别投10000次的结果，以频率代表投骰子的概率。再给出一组可能的结果，问是否有一种投骰子的方式，使得投10000次后的结果刚好与给出的相符合，而且要求每个骰子至少投一次。很显然，我们很快的可以把问题转化为求一个非齐次线性方程组是否有正整数解的问题。设第一个骰子为a，对应三个面的概率分别为a1、a2、a3，第二个骰子为b、第三个为c。设一组可能的解向量为[x，y，z]，最后的结果为[h1，h2，h3]。那么就相当于看方程组ai*x+bi*y+ci*z=hi、x+y+z=10000是否有正整数解。

        由于个人的脑残，我一开是并没有考虑xyz之和为10000，也没有考虑当矩阵不满秩的时候的解的情况(认为既然是无穷解，肯定有符合要求的整数解)，所以赛场上疯狂的WA。我都不想说自己了……

        场下，我认识到应该是四个方程求解，根据解的情况分情况讨论：

                1、无解。无解的情况自然不存在，输出NO。

                2、存在唯一解。只需要判断这个唯一解是否正整数即可。

                3、存在无穷解。这里又要分几种小的情况：

                            a、存在一个自由变元。这意味着矩阵的秩为2。那么，我们只需要枚举这个自由变元的值(1~9998)，然后对应算出另外两个变量的值，看是否满足条件即可

                            b、存在两个自由变元。由于变量至少满足方程x+y+z=10000，而又有两个自由变量，故此情况一定存在正整数解

        就是这样没有更多的情况了。然后对应的用高斯消元即可。高斯消元的具体步骤在代码中会有解释，代码如下：

#include <bits/stdc++.h>
#define eps 1e-9

using namespace std;

double a[4][4];
int x[3];
bool free_x[3];


int Gauss(int equ,int var)
{
	int i,j,k,max_r,col;
	int free_x_num,free_index;
	double temp,lcm,t;

	for(int i = 0;i < var;i++)			//初始化
	{
		x[i] = 0;				//开始默认全为自由变元
		free_x[i] = true;
	}

	col = 0;
	for(k = 0;k<equ && col<var;k++,col++)
	{
		max_r = k;
		for(i = k+1;i < equ;i++)		//找到对应变量的最大的系数，减小相除时产生的误差
			if(a[i][col]-a[max_r][col]>eps) max_r = i;
		if(max_r!=k)
			for(j = k;j < var+1;j++) swap(a[k][j],a[max_r][j]);	//初等行变换
		if(fabs(a[k][col])<eps) {k--; continue;}
		for(i = k+1;i < equ;i++)
			if(a[i][col]!=0)
			{
				t = a[i][col]/a[k][col];
				for(j = col;j < var+1;j++)
					a[i][j] = a[i][j]-a[k][j]*t;		//对应列的系数消为0，和线性代数的解法相同
			}
	}

	for(i = k;i < equ;i++)
		if(fabs(a[i][col])>eps) return -1;			//如果出现系数为0，但方程右侧不为0，则说明无解
	if(k < var)
	{
		for(i = k-1;i>=0;i--)				//已经是上三角矩阵，开始尝试求解
		{
			free_x_num = 0;
			for(j = 0;j < var;j++)			//寻找系数不为0的变量
				if(fabs(a[i][j])>eps && free_x[j]) free_x_num++,free_index = j;
			if(free_x_num > 1) continue;		//变量数目大于1，说明有自由变元
			temp = a[i][var];
			for(j = 0;j < var;j++)
				if(fabs(a[i][j])>eps && j!= free_index) temp -= a[i][j]*x[j];	//进行移项
			x[free_index] = temp/a[i][free_index];		//求出非自由变元的解
			if (fabs(x[free_index])<eps) return -3;		//负数解，不符合要求
			if (fabs(x[free_index]-int(x[free_index]))>eps) return -2;	//浮点数解，不符合要求
			free_x[free_index] = 0;			//标记非自由变元
		}
		return (var-k);					//返回自由变元的个数
	}
	for(i = var-1;i>=0;i--)
	{
		temp = a[i][var];
		for(j = i+1;j<var;j++)
			if(fabs(a[i][j])>eps) temp -= a[i][j]*x[j];
		x[i] = temp/a[i][i];				//迭代继续求解
	}
	return 0;						//出现唯一解
}


int main()
{
    while (1)
    {
        for(int i=0;i<3;i++)
        {
            int s=0;
            for(int j=0;j<3;j++)
            {
                scanf("%lf",&a[j][i]);
                s+=a[j][i];
                a[j][i]/=10000;
            }
            if (s==0) return 0;
        }
        for(int i=0;i<3;i++)
        {
            scanf("%lf",&a[i][3]);
            a[3][i]=1;
        }
        a[3][3]=10000;
        int num=Gauss(4,3),flag=0;
        if (num==0)
        	for(int i=0;i<3;i++)
        		if (x[i]<=0) flag=1;
        if (num==1)					//自由变元数目为1时
        {
        	double x,y,z,s;
        	x=a[1][0],y=a[1][1];
        	z=a[1][2],s=a[1][3];
        	flag=1;
        	if (fabs(x)<eps)			//暴力枚举自由变元，判断是否存在正整数可行解
        		for(int i=1;i<int(s/y+0.5);i++)
        		{
        			double k=(s-i*y)/z;
        			if (fabs(k-int(k+0.5))<eps&&k+i<10000) {flag=0;break;}
				}
		if (fabs(y)<eps)
        		for(int i=1;i<int(s/x+0.5);i++)
        		{
        			double k=(s-i*x)/z;
        			if (fabs(k-int(k+0.5))<eps&&k+i<10000) {flag=0;break;}
				}
		if (fabs(z)<eps)
        		for(int i=1;i<int(s/y+0.5);i++)
        		{
        			double k=(s-i*y)/x;
        			if (fabs(k-int(k)+0.5)<eps&&k+i<10000) {flag=0;break;}
				}
		}
        if (num<0||flag) puts("NO"); else puts("YES");	
    }
    return 0;
}