matlab 및 C++ 매트릭스 제품 : 동일한 입력, 다른 출력

Question

x86 컴퓨터의 우분투 12.04 및 g ++ 4.6.3에서 Matlab 2010을 사용하고 있습니다. 내가 생산과 입력 어떻게 입니다 :

#include <Src/Tools/Math/Matrix_nxn.h> 
    #include <iostream> 
    using namespace std; 
    int main() 
    { 
     Matrix_nxn<double,4> A1,A2,Tb,aa; 
     A1[0][0] = 0.99958087959447828; A1[0][1] = 1.7725781974830023e-18;A1[0][2] = 0.028949354900049871; A1[0][3] = 0; 
     A1[1][0] = -0.028949354900049871; A1[1][1] = 6.1204654815537932e-17;A1[1][2] = 0.99958087959447828; A1[1][3] = 0; 
     A1[2][0] = 0,   A1[2][1] = -1;   A1[2][2] = 6.1230317691118863e-17;A1[2][3] = 0.21129000000000001; 
     A1[3][0] = 0,   A1[3][1] = 0;   A1[3][2] = 0;    A1[3][3] = 1; 

     A2[0][0] = 0.90634806393366396; A2[0][1] = -0.42253187690835708;A2[0][2] = 0;A2[0][3] = 0; 
     A2[1][0] = 0.42253187690835708; A2[1][1] = 0.90634806393366396; A2[1][2] = 0;A2[1][3] = 0; 
     A2[2][0] = 0;   A2[2][1] = 0;   A2[2][2] = 1;A2[2][3] = 0; 
     A2[3][0] = 0;   A2[3][1] = 0;   A2[3][2] = 0;A2[3][3] = 1; 

     Tb[0][0] = 0.99956387949834924; Tb[0][1] = -0.00016363183229951183; Tb[0][2] = -0.029530052943282908; Tb[0][3] = 0; 
     Tb[1][0] = 0;   Tb[1][1] = 0.99998464792303143; Tb[1][2] = -0.0055411116439683869;Tb[1][3] = 0; 
     Tb[2][0] = 0.029530506297888514;Tb[2][1] = 0.0055386950515785164; Tb[2][2] = 0.99954853411673616; Tb[2][3] = 0; 
     Tb[3][0] = 0;   Tb[3][1] = 0;   Tb[3][2] = 0;    Tb[3][3] = 1; 


    aa = Tb*A1*A2; 
    cout.precision(25); 
    cout <<aa[0][0]<<' '<<aa[0][1]<<' '<<aa[0][2]<<' '<<aa[0][3]<<endl 
     <<aa[1][0]<<' '<<aa[1][1]<<' '<<aa[1][2]<<' '<<aa[1][3]<<endl 
     <<aa[2][0]<<' '<<aa[2][1]<<' '<<aa[2][2]<<' '<<aa[2][3]<<endl 
     <<aa[3][0]<<' '<<aa[3][1]<<' '<<aa[3][2]<<' '<<aa[3][3]<<endl; 
} 

을하고이 operator*의 정의입니다 :

이

Matrix_nxn<T, N> res; 
size_t i, j, k; 
for (i = 0; i < N; ++i) 
{ 
    for (j = 0; j < N; ++j) 
    { 
    for (k = 0; k < N; ++k) 
    { 
     res[i][j] += m1[i][k] * m2[k][j]; 
    } 
    if (MVTools::isNearInf(res[i][j])) 
    { 
     if (MVTools::isNearPosInf(res[i][j])) 
     throw MVException(MVException::PosInfValue); 
     else 
     throw MVException(MVException::NegInfValue); 
    } 
    } 
} 
return res; 

이상한 것은 내가 MATLAB 및 I 안에 같은 값으로 같은 행렬을 만드는 것입니다 다른 결과를 얻는다. 다음은 Matlab 코드입니다.

Tb = [0.99956387949834924,-0.00016363183229951183,-0.029530052943282908,0;0,0.99998464792303143,-0.0055411116439683869,0;0.029530506297888514,0.0055386950515785164,0.99954853411673616,0;0,0,0,1]; 
A1 = [0.99958087959447828,1.7725781974830023e-18,0.028949354900049871,0;-0.028949354900049871,6.1204654815537932e-17,0.99958087959447828,0;0,-1,6.1230317691118863e-17,0.21129000000000001;0,0,0,1]; 
A2 = [0.90634806393366396,-0.42253187690835708,0,0;0.42253187690835708,0.90634806393366396,0,0;0,0,1,0;0,0,0,1]; 
aa = Tb*A1*A2; 
aa - aaa 

ans = 

    1.0e-16 * 

       0 -0.555111512312578     0     0 
       0     0     0     0 
       0     0     0     0 
       0     0     0     0 

aaa는 C++ 구현의 출력입니다. 오류가 너무 작다는 것을 알고 있지만 문제의 원인을 알고 싶습니다! 많은 코드를 디버깅하고 싶고 좋은 디버깅을 위해서는 0 차이가 필요합니다.

Answer 1

C++ 코드에서 25 자리의 정밀도를 기대합니다. double 유형을 사용하는 것은 거의 불가능합니다. long double을 사용하면 더 나은 정밀도를 얻을 수 있지만 25 자리와 같은 음흉은 아닙니다.

는

참조 : What is the precision of long double in C++?

Answer 2

다른 값에 대한 이유는 (그러나 그것은 수 있습니다 미미)는 MATLAB에서 사용하는 알고리즘과 그렇지 같은 을 점이다.

알고리즘은 간단합니다. O(N^3) 행렬 곱셈. 효율적으로 계산할 수있는 작은 크기의 행렬에 대한 특별한 알고리즘과 O(N^3)보다 나은 점근 적 동작을 갖는 복잡한 알고리즘이 있습니다. 소위 작은 오류의

• Laderman's 3x3 matrix multiplication with only 23 multiplications, is it worth it?
• http://en.wikipedia.org/wiki/Strassen_algorithm