高斯消元

      话说今天模拟赛P1（就是无向图的生成树个数一文中的例题囧||）是一个极其强大的数学问题，要用到矩阵这种让人一看就头晕的东西，而saltless在看完黑客帝国那个名字叫做”矩阵革命”的纠结电影后，一直对矩阵、行列式等东西很是恶心。谁料到，P2竟然又是一个用到行列式的高斯消元……恶心啊～

      高斯消元是干什么的呢？对于n元一次方程组，我们可以用手工的方法进行逐个消元。一想那一堆方程，头皮都麻了，编程实现就更找不到头绪了。高斯消元就是用短短的一小段代码来解n元一次方程组。

      高斯消元怎么做呢？听着名字很强悍，比如说，名字很强悍的欧拉函数就是一个极其强大的东西…其实，高斯消元就是像小学讲的高斯在小石板上写下的5050一样，是一个很简单的东西。

      回忆下小学和初中时解二元一次不等式的过程。我们先找到x的两个系数的最小公倍数，然后将x的系数扩大到一样，这样就能够消掉x，求得y的值，再将y代回方程就可以解出x。

      其实，高斯消元就是这样的一个过程。

      我们都知道行列式的代数性质（注意，第n+1列存储本方程等号右侧的系数）。这样，把第一行的要消去的元的系数和下面几行的相应的元的系数通过放缩化成一致（由于是实数，只放缩其中一个系数即可），进行减法即可消去该元。然后依次用第i行消去[i+1,n]的xi，直到只剩两项，即xn的系数和这时的等式右边的常数kn。只要用kn/xn就可以求出xn的近似解（在一定的数据范围内，double或者extended的精确度足够，这一点在下文会涉及），然后依次回代，即可求出所有x值。

      听着很玄乎，配合代码看就会很简单。

      另外说一些关于精度的问题。实数类型不可避免的涉及精度的问题。显然，当n很大时，不断的进行除法运算会将近似更加近似，影响效果。这时候我们可以用整数进行运算，也就是找到两个数的最小公倍数进行减法就可以避免除法的近似。但是，整型会出现被0除的错误，需要特殊处理，并且当要消的元两个系数是两个大质数的话，行列式中的元素就会变得很大。如果再进行化简，就是又一重时间和空间的浪费，并且编程复杂度大大增高。所以，在一定的数据范围内，用实数完全可以实现高斯消元。

      下面是一道例题。

 

【问题描述】
    终于逃出了夜爸爸的魔掌，小夜来到了一个巨大无比的超市。超市里琳琅满目的商品勾起了小夜
的购买欲， 于是小夜推来n辆大型购物车••开始往里面塞东西！@#￥……塞了一车又一车……因为太
呆了，小夜喜欢重复地买一些东西，而且每辆购物车里都会有n件物品。比如她在推第一辆车的时候会
买熊宝宝、小saber、机器人……在推第2辆车的时候她还是会买熊宝宝、小saber、机器人，而且买
这些东西的顺序一样，只是数目不同罢了……结帐的时候，售货员告诉她每辆购物车里的商品总价值
x元，累加起来后小夜发现这是个天文数字……她开始怀疑售货员有没有坑她。因为每辆车里商品太多，
小夜不能直接看到商品的单价，她只知道第i辆购物车里第ai件商品的数目是多少。但她根本算不出
那么多物品的单价，于是她想到了你……
    任务：小夜给了你一张表，上面写着每辆车里各个商品的数目以及商品总价，请你帮小夜算出每
件商品的单价是多少。
【输入格式】
    输入文件中第一行有一个数字n代购物车数
    接下来有n行，每行n+1个整数。前n个整数aij  代表第i辆辆购物车里第j件商品的数目，第n+1
个整数xi代表该辆购物车的总价值。
【输出格式】
输出文件一共有n行，每行一个整数，代表第i件物品vi的单价。
【输入样例】
2
1 1 3
1 2 5
【输出样例】
1
2
【数据范围】
对于20%的数据： 0<n<=10
对于100%的数据：  0<n<=500  0<vi<32767  0<xi<32767   0<aij<101
对于10%的数据限1s，对于另外20%的数据时限2s

 

 

      首先注意到，数据范围比较小。用实数进行高斯消元解方程组即可。

 

参考代码（递归）：

 

program shopping;
  var
    a:array[1..501,1..501]of double;
    i,j,n:integer;
    k:array[1..500]of double;
  procedure kill(x:integer);  //消元
    var
      i,j:integer;
      zoom,tot:double;
    begin
      if x=n then
        begin
          k[n]:=a[n,n+1]/a[n,n];
          exit;
        end;
      for i:=x+1 to n do
        begin
          zoom:=a[i,x]/a[x,x];  //求放缩的比例
          for j:=x+1 to n+1 do
          a[i,j]:=a[i,j]-a[x,j]*zoom;  //消元
        end;
      kill(x+1);
      tot:=0;
      for i:=x+1 to n do
        tot:=tot+k[i]*a[x,i];  //将值代回方程，求得上一个未知数的值
      k[x]:=(a[x,n+1]-tot)/a[x,x];  
    end;
  begin
    readln(n);
    for i:=1 to n do
      for j:=1 to n+1 do
        read(a[i,j]);
    kill(1);
    for i:=1 to n do
      writeln(k[i]:0:0);
  end.

 

参考代码（非递归）：

 

program shopping;
  var
    a:array[1..501,1..501]of double;
    i,j,n,k:integer;
    zoom,sum:double;
  begin
    readln(n);
    for i:=1 to n do
      for j:=1 to n+1 do
        read(a[i,j]);
    for k:=1 to n-1 do
      for i:=k+1 to n do
        begin
          zoom:=a[i,k]/a[k,k];
          for j:=k+1 to n+1 do
            a[i,j]:=a[i,j]-a[k,j]*zoom;
        end;
    for i:=n downto 1 do
      begin
        sum:=0;
        for j:=n downto i+1 do
          sum:=sum+a[N+1,j]*a[i,j];
        a[n+1,i]:=(a[i,n+1]-sum)/a[i,i];
      end;
    for i:=1 to n do
      writeln(a[n+1,i]:0:0);
  end.

 

 

本文地址：http://www.cnblogs.com/saltless/archive/2010/11/07/1870928.html

（saltless原创，转载请注明出处）