这年头，还有人用C++的吗？

最近在时间的细屑里面挤出病毒般大小的时间来看GF调不出来的程序，简直是郁闷死人了。又要看好多年没有真正用过的C++，真是一种头痛欲裂的感觉，在一次感受到C++的疯狂之处……（这个程序其实之前确实是在用C++写的，后来因为某种原因，放弃了C++当中的++部分，用了几乎纯粹C的功能来编写。准确一点讲，标题应该是“这年头，还有人用C的吗？”）

大家先来看两张对比的图，着一张是正确的输出结果：


再来看一张错误的输出结果：（这一个是我需要调的程序）



然后就开始要求我GF赶紧对公式，对比正确输出所用程序的代码，看看两者到底有什么地方不一样。由于两个程序所用的坐标系不一样，所以对比起来比较困难，最后费了好大的劲，把坐标系修改成完全一样的之后，发现问题依旧。这个时候另外一个现象更加困惑了我：在修改错误程序的坐标系之后，发现图形的缩小程度明显改善了，但是仍然是严重收敛的。这个现象着实打击了我一下，因为我一开始只是以为是误差量的积累造成的。难道说就是调整了一下几个公式的位置以及一些正负号就会造成误差累积量的不同？还真是前所未见的情况！
这里整个调试过程是非常痛苦的，因为总共大约有几千万次的循环，计算量之大不可想象。而输出正确和错误的程序之间，在大约每60次循环才会产生大约1x10 ^-7细微差别。难为我GF在那里单步调试了几十万次的循环，真是可怜啊。到底是哪里出了问题呢？多亏我用VC7来运行了一遍，才发现问题所在。（之前都是我在抽空Review代码，GF在1000km+的距离以外用VC6在调试。）下面这段代码就是存在错误的地方：

void  MultiMatr( int  s,  int  n,  int  t,  double   * arr1,  double   * arr2,  double   * result)
{
    double *tmp;

    for (int i = 0; i < s; i++)
    {
        for (int j = 0; j < t; j++)
        {
            tmp = result + t * i + j;
            *tmp = 0;
            for (int k = 0; k < n; k++)
            {
                *tmp += ((*(arr1 + n * i + k)) * (*(arr2 + t * k + j))); 
            }
        }
    }
}

/**/ /***************************** update ***************************/
void  UpdateQ( double  WnbbA[ 3 ])
{
    int        i;
    double    deltaSita0, angSin, angCos;
    double    deltaSita[3], qChange[3], qTrans[4][4];

    deltaSita0 = 0;
    
    for (i = 0; i < 3; i++)
    {
        deltaSita[i] = WnbbA[i] * quaDel;
        deltaSita0 += deltaSita[i] * deltaSita[i];
    }
    deltaSita0 = sqrt(deltaSita0);

    angCos = cos(deltaSita0 / 2);
    if (deltaSita0 == 0)
        angSin = 0.5;
    else
        angSin = sin(deltaSita0 / 2.0)/deltaSita0;

    qTrans[0][0] = qTrans[1][1] = qTrans[2][2] = qTrans[3][3] = angCos;
    qTrans[0][1] = qTrans[3][2] = -(angSin * deltaSita[0]);
    qTrans[1][0] = qTrans[2][3] = -qTrans[0][1];
    qTrans[0][2] = qTrans[1][3] = -(angSin * deltaSita[1]);
    qTrans[2][0] = qTrans[3][1] = -qTrans[0][2];
    qTrans[0][3] = qTrans[2][1] = -(angSin * deltaSita[2]);
    qTrans[3][0] = qTrans[1][2] = -qTrans[0][3];

    MultiMatr(4, 4, 1, (double*)qTrans, (double*)q, (double*)qChange);

    for (i = 0; i < 4; i++)
        q[i] = qChange[i];

}

大家看出来那里错了吗？真是稀松平常的错误啊！害死人了。