(((silence)))

// Author: IF

// License: BSD

double determinant( NA::MatDP const& m )
{
    NA_ASSERT(m.rows()==m.cols());
    NA::MatDP mat(m);
    double    d = 1;
    int       i,j,k,n = m.rows();

    for( i=0; i<n; ++i ) {
        double largest=fabs(mat[i][i]);
        int    largestAtRow=i;
        // 选取主元
        for ( j=i+1; j<n; ++j ) {
            if ( fabs(mat[j][i]) > largest ) {
                largest = fabs(mat[j][i]);
                largestAtRow = j;
            }
        }
        if ( largest < 1e-20 ) // 0
            return 0;
        // 行交换
        if ( largestAtRow != i ) {
            mat.swapRow ( i, largestAtRow );
            d = -d;
        }
        // 消去
        for ( j=i+1; j<n; ++j ) {
            double factor = mat[j][i]/mat[i][i];
            mat[j][i]=0;
            for ( k = i+1; k<n; ++k ) {
                mat[j][k] -= mat[i][k]*factor;
            }
        }
    }
    // 计算行列式
    for ( i=0; i<n; ++i ) {
        d *= mat[i][i];
    }

    return d;
}

// 测试：

int main()
{
    double pm[]= {
        1, 5, 5,
        3, 2, 2,
        5, 5, 8
    };
    NA::MatDP m(pm, 3,3);
    cout<<determinant(m)<<endl;
    return 0;
}