求一数学问题解决办法：矩阵的特征向量

Question

现在有个矩阵计算的需求，需要求出一实数范围上的n阶方阵的特征向量，其中n为常数（即每次输入的方阵阶数都相同），希望下面几个问题可以得到回答：
1. 是否有与此矩阵计算相关的php内部函数？
2. 是否有已经知道的前人写过的相关的函数？
3. 如果上面两个都没有，那么是否可以通过其他语言的程序（比如一个C++程序）来计算，如果是，用php怎么粘接？
4. 如果有人愿意写一个，再感激不过了（额外加100分）

Answer 1

这是百度百科里的特征向量词条
http://baike.baidu.com/view/475996.html我也是昨天才翻出线性代数书复习了概念，上面只介绍了一种解法，需要求解齐次线性方程组定义：设A是数域F上的一个n阶方阵，如果存在 λ∈F，使 AX=λX 有非零解：X=(ξ1,ξ2,...,ξn)'∈Fn，则称λ为矩阵A的一个特征值，而称X为属于特征值的一个特征向量。

Answer 2

 现将2002年用构造了一个矩阵类，其中有编写了一方法，通过Jocobbi迭代计算特征根与特征向量，建议楼主将其转换为php函数，应该不难。
矩阵v为特征向量矩阵， pArray为为特征值
------------------------------------------------
void Matrix :: Jacobbi ( Matrix& v , double* pArray )
{
  int p , q , j , ind ,n;
  double dsqr , d1 , d2 , thr , dv1 , dv2 , dv3 , dmu , dga , st , ct ;
  double eps = 0.00000001;
  int* iZ;  //add 2002.8.27  Matrix CA ( *this ) ;
  n = Rownum ;  //add 2002.8.27
  try {
    iZ = new int[n] ;
  }
  catch (std::bad_alloc) {  // ENTER THIS BLOCK ONLY IF bad_alloc IS THROWN.
      // YOU COULD REQUEST OTHER ACTIONS BEFORE TERMINATING
    Application -> MessageBox("不能在堆中分配内存,再见...","异常", MB_OK);
    exit(-1);
  }
  //add 2002.8.27  for ( p = 0 ; p < n ; p ++ )
    for ( q = 0 ; q < n ; q ++ )
      v.dmat[p][q] = ( p == q ) ? 1.0 : 0 ;  dsqr = 0 ;
  for ( p = 1 ; p < n ; p ++ )
    for ( q = 0 ; q < p ; q ++ )
      dsqr += 2 * CA.dmat[p][q]* CA.dmat[p][q] ;
  d1 = sqrt ( dsqr ) ;
  d2 = eps / n * d1 ;
  thr = d1 ;
  ind = 0 ;
  do {
    thr = thr / n ;
    while ( !ind ) {
    for ( q = 1 ; q < n ; q ++ )
      for ( p = 0 ; p < q ; p ++ )
if ( fabs ( CA.dmat[p][q] ) >= thr ) {
  ind = 1 ;
  dv1 = CA.dmat[p][p] ;
  dv2 = CA.dmat[p][q] ;
  dv3 = CA.dmat[q][q] ;
  dmu = 0.5 * ( dv1 - dv3 ) ;
          double dls = sqrt ( dv2 * dv2 + dmu * dmu ) ;
  if ( fabs ( dmu ) < 0.00000000001 )  dga = -1 ;
          //if ( dmu == 0.0 )  dga = -1.0 ;
  else dga = ( dmu < 0 ) ? ( dv2 / dls ): ( - dv2 / dls );
          st = dga / sqrt ( 2 * ( 1 + sqrt ( 1 - dga * dga ) ) ) ;
  ct = sqrt ( 1 - st * st ) ;
          for ( register int l = 0 ; l < n ; l ++ ) {
    dsqr = CA.dmat[l][p] * ct - CA.dmat [l][q] * st ;
    CA.dmat [l][q] = CA.dmat[l][p] * st + CA.dmat[l][q] * ct ;
    CA.dmat [l][p] = dsqr ;
    dsqr = v.dmat[l][p] * ct - v.dmat[l][q]*st ;
    v.dmat[l][q] = v.dmat[l][p] * st + v.dmat[l][q] * ct ;
    v.dmat[l][p] = dsqr ;
  }
  for ( register int l = 0 ; l < n ; l ++ ) {
    CA.dmat[p][l] = CA.dmat[l][p] ;
    CA.dmat[q][l] = CA.dmat[l][q] ;
  }
  CA.dmat [p][p] = dv1 * ct * ct + dv3 * st * st - 2 * dv2 * st * ct;
  CA.dmat [q][q] = dv1 * st * st + dv3 * ct * ct + 2 * dv2 * st * ct;
          CA.dmat [p][q] = CA.dmat[q][p] = 0.0 ;
        }
if ( ind ) ind = 0 ;
else break ;
      }
  } while ( thr > d2 ) ;
  for ( register int l = 0 ; l < n ; l ++ ) {
    pArray[l] = CA.dmat[l][l];
    iZ [l] = l;
  }
  //2002.8.27 add
  double dTemp ;
  int i , k;  for ( i = 0 ; i < n ; i++ ){
    //dmax = pArray[i];
    for ( j = i + 1 ; j < n ; j++ ){
      if ( pArray[i] < pArray[j]){
        dTemp = pArray[i];
        pArray[i] = pArray[j];
        pArray[j] = dTemp;
        k = iZ[i] ;
        iZ[i] = iZ[j];
        iZ[j] = k;
      }
    }
  }
  CA = v;  for ( j = 0 ; j < n ; j++ )
    for ( i = 0 ; i < n ; i++ )
      v.dmat [i][j] = CA.dmat[i][iZ[j]];  delete[] iZ;
  //2002.8.27 end add
}