那天志旭看到上海赛区的题目挂到spoj上面了. 我也突然来了兴致, 想看看这道题目是不是和我想的一样. 说实话,这道题目其实就是我们校赛预赛的那道三维照相机的题的简化版本+最大独立集.

比赛现场我还在想我们学校的那两个队伍应该有最少有一个会搞这道题. 没有想到他们都怕40个点的最大独立集会超时...不敢去做这道题目. 真是faint. 碰上个原题, 居然不敢做.

很明显这道题目的数据很难出的很强. 才40个点, 明显就是让你暴力过去的.

嗯...不废话了. 说说题目解法.

几何部分的核心就是求出三维三角形的遮挡关系.

这一部分可以这样求: 对于每一个三角形, 从原点开始对三个点分别做三条射线. 这样就包含出了一个无限的空间. 那么怎么判断这些无限空间是不是有相交呢?

一对空间的3 * 3 == 9对"面"如果有相交的话, 那么这对空间就一定相交了. 否则只剩下包含的情况.

包含的情况只要在某个空间取出一个点, 看看是不是在另一个空间之内就好了.

 矩阵真的好优美，好多在数学中那些乘乘加加的操作，其本质都是矩阵。这些操作，往往都是把事物进行分解，对每一个基量进行操作，最后再把这些操作整合起来。

向量旋转，这个几何中的基本操作，可以相当优美的用矩阵来操作。究其原因，也是因为向量是多维的。用矩阵来操作，正是把向量分解开来，分别旋转，最后再进行整合。

看下面的矩阵

这个是二维向量的旋转矩阵，它可以将一个向量逆时针旋转一个角度。

将其变形，变会得到二维向量的顺时针旋转的形式：

这里要注意一种特殊的情况，就是当角度为90度的时候，sin和cos的结果只有1, 0, -1三种可能性。所以这个矩阵可以改写成特殊的形式，其意义在于用这样的旋转操作，不会产生精度问题。sin和cos的运算的精度是比较低的，能少用则尽量少用。在计算几何中的向量旋转操作，大部分都可以通过变形，只用到旋转90度，从而避免精度问题。

再来说一点就是，这些旋转矩阵都有一些特点，最明显的莫过于他们的行列式的值都是1。所以在验证正确性的时候，可以用这个操作。也正因为有这一个性质，所以向量才只会进行旋转，不会进行缩放。

三维向量的旋转矩阵：

上面三个公式的旋转方向可以看成按右手定则。

来看第三个公式，这个公式是围绕Z轴，把X轴往Y轴方向移动。我们把拇指向上（表示Z轴），手指指向X轴，然后手指自然弯曲方向便是旋转方向。

其它两个公式也是类似。

其实第三个公式，去掉Z轴，就是开头讲的逆时针旋转的公式。

在其它方向上的旋转都可以由这三个矩阵组合而来。

围绕轴u = (ux, uy, uz)来旋转的矩阵代码如下：

double a[3][3] = {
        {SQR(ux) + (1 - SQR(ux)) * c, ux * uy * (1 - c) - uz * s, ux * uz * (1 - c) + uy * s},
        {ux * uy * (1 - c) + uz * s, SQR(uy) + (1 - SQR(uy)) * c, uy * uz * (1 - c) - ux * s},
        {ux * uz * (1 - c) - uy * s, uy * uz * (1 - c) + ux * s, SQR(uz) + (1 - SQR(uz)) * c}
    };

要注意这个ux * ux + uy * uy + uz * uz = 1

有一个比较有意思的问题就是，知道旋转矩阵之后，怎么来确定向量u呢？

我们做如下变形之后，可以得到：

Ru = Iu –> (R – I)u = 0

也就是说我们要找到一个非0的向量，使得他和一个矩阵乘起来得到一个空矩阵。这个问题可以用高斯消元来解决。实际上我们就是要解一个线性方程组。

 题目连接：http://www.topcoder.com/stat?c=problem_statement&pm=10516&rd=13900

昨天做SRM再次失败，只做出了第一道题。每次做500分都差一点。郁闷的要死，昨天又是死活不过样例，今天醒来想一想，有一个细节方面还没有想清楚就急于去写了。导致算法上有很大的纰漏。

实际上这题算法还是挺容易想的，不过前提是要对矩阵乘法比较了解。实际上这道题也是代表了矩阵乘法的一种类型。我们比较常接触的是那种用乘法和加法的矩阵。实际如下替换之后，矩阵的各种性质也是保持的：

乘法-->加法

加法-->取最大值

这种矩阵，我们的乘法操作是如下的：

 //假设我们是要把Matix A * Matix B的结果存到Matix ans中
 for (i = 0 to n  -1 )
      for (j = 0 to n – 1)
          ans[i,j] = minValue;
          for (k = 0 to n – 1)
                 ans[i][j] = max(ans[i][j], a[i][k] + b[k][j])
从这种操作的意义上来看，就是用来取最长路之类的。跟本题是一样的。

这个题目特殊的地方在于，它有两个要求，一种是必须走K步，中间不可以停下，另一种就是中间可以停。

比赛时我没有看清楚这个条件，以为都是可以停的，所以直接按最原始的做法来做的，等我反应过来看错题的时候，又太急了，没有想清楚就乱去写了，导致算法写反了。

算法如下:

首先我们构造原始矩阵A，A就是由输入的那三个参数来的，但是要加一种情况就是如果两点不可达的时候，要把值设置成无穷小。而不是题目里的0.

然后我们求A = A^stepsPerSecond。这样便可以求出走这些步之后的情况。

再下一步我们不可以直接求A^timeLimit，为什么呢？这是因为我们中间是可以停下的。为了处理这种情况，我们可以把A[1,1]设置为0。这样我们便可以求得最终的答案了。不过这里要注意A[1,1]有可能本来是大于0的，这个时候就不用管它了。不要把这里也设成0。

这种矩阵还有一点和普通的不同在于它的单位矩阵并不是对角线全是1。而是对角线全是0，其它的全是无穷小。这样的矩阵和另一个矩阵相乘之后才可以保持另一个矩阵的不变。

上面如果没有看明白，就直接看下面的代码吧，代码还是比较清楚的：

public class AntOnGraph
{
    static int n;
    public string maximumBonus(string[] p0, string[] p1, string[] p2, int stepsPerSecond, int timeLimit)
    {
        n = p0.Length;
        long[,] st = new long[n, n];
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                int k = (int)(p0[i][j] - '0') * 100 + (int)(p1[i][j] - '0') * 10 + (int)(p2[i][j] - '0');
                k -= 500;
                if (k == -500) {
                    st[i,j] = long.MinValue;
                } else {
                    st[i,j] = k;
                }
            }
        }
        st = mypow(st, stepsPerSecond);
        if (st[1,1] < 0) {
            st[1,1] = 0;
        }
        st = mypow(st, timeLimit);
        if (st[0,1] == long.MinValue) {
            return "IMPOSSIBLE";
        } else {
            return st[0, 1].ToString() ;
        }
    }
    public long[,] mypow(long[,] a, int b)
    {
        long[,]ans = new long[n, n];
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                if (i == j) {
                    ans[i,j] = 0;
                } else {
                    ans[i,j] = long.MinValue;
                }
            }
        }
        for (; b != 0; b >>= 1) {
            if (b % 2 != 0) {
                ans = mul(ans, a);
            }
            a = mul(a, a);
        }
        return ans;
    }
    public long[,] mul(long[,] a, long[,] b) {
        long[,] c = new long[n, n];
        for (long i = 0; i < n; ++i) {
            for (long j = 0; j < n; ++j) {
                c[i,j] = long.MinValue;
                for (int k = 0; k < n; ++k) {
                    if (a[i,k] > long.MinValue && b[k,j] > long.MinValue) {
                        if (a[i,k] + b[k,j] > c[i,j]) {
                            c[i,j] = a[i,k] + b[k,j];
                        }
                    }
                }
            }
        }
        return c;
    }
}