Простой 3x3 матричный обратный код (C++)

Question

Простой 3x3 матричный обратный код (C++)

Попробуйте это:

.box {
  outline: 2px solid #ddd;
  margin-top: -2px;
  border-top: 10px solid #3F9BD0;
  min-width:100px;
  min-height:100px;
  float:left;
}

<div class="box"></div>

36

matrix c++ math matrix-inverse

задан genpfault 15 February 2012 в 23:12

6 ответов

При всем уважении к нашему неизвестному (yahoo) плакату, я смотрю на такой код и просто умираю немного внутри. Алфавитный суп безумно сложно отладить. Единственная опечатка где-нибудь там действительно может испортить вам весь день. К сожалению, в этом конкретном примере отсутствуют переменные с подчеркиванием. Намного веселее, когда у нас есть a_b-c_d * e_f-g_h. Особенно при использовании шрифта, где _ и - имеют одинаковую длину в пикселях.

Принимая предложение Сувеша Пратапы, я отмечаю:

Given 3x3 matrix:
       y0x0  y0x1  y0x2
       y1x0  y1x1  y1x2
       y2x0  y2x1  y2x2
Declared as double matrix [/*Y=*/3] [/*X=*/3];

(A) Когда мы берем минор из массива 3x3, мы получаем 4 представляющих интерес значения . Нижний индекс X / Y всегда равен 0 или 1. Более высокий индекс X / Y всегда равен 1 или 2. Всегда! Следовательно:

double determinantOfMinor( int          theRowHeightY,
                           int          theColumnWidthX,
                           const double theMatrix [/*Y=*/3] [/*X=*/3] )
{
  int x1 = theColumnWidthX == 0 ? 1 : 0;  /* always either 0 or 1 */
  int x2 = theColumnWidthX == 2 ? 1 : 2;  /* always either 1 or 2 */
  int y1 = theRowHeightY   == 0 ? 1 : 0;  /* always either 0 or 1 */
  int y2 = theRowHeightY   == 2 ? 1 : 2;  /* always either 1 or 2 */

  return ( theMatrix [y1] [x1]  *  theMatrix [y2] [x2] )
      -  ( theMatrix [y1] [x2]  *  theMatrix [y2] [x1] );
}

(B) Определитель теперь: (Обратите внимание на знак минус!)

double determinant( const double theMatrix [/*Y=*/3] [/*X=*/3] )
{
  return ( theMatrix [0] [0]  *  determinantOfMinor( 0, 0, theMatrix ) )
      -  ( theMatrix [0] [1]  *  determinantOfMinor( 0, 1, theMatrix ) )
      +  ( theMatrix [0] [2]  *  determinantOfMinor( 0, 2, theMatrix ) );
}

(C) И обратное теперь:

bool inverse( const double theMatrix [/*Y=*/3] [/*X=*/3],
                    double theOutput [/*Y=*/3] [/*X=*/3] )
{
  double det = determinant( theMatrix );

    /* Arbitrary for now.  This should be something nicer... */
  if ( ABS(det) < 1e-2 )
  {
    memset( theOutput, 0, sizeof theOutput );
    return false;
  }

  double oneOverDeterminant = 1.0 / det;

  for (   int y = 0;  y < 3;  y ++ )
    for ( int x = 0;  x < 3;  x ++   )
    {
        /* Rule is inverse = 1/det * minor of the TRANSPOSE matrix.  *
         * Note (y,x) becomes (x,y) INTENTIONALLY here!              */
      theOutput [y] [x]
        = determinantOfMinor( x, y, theMatrix ) * oneOverDeterminant;

        /* (y0,x1)  (y1,x0)  (y1,x2)  and (y2,x1)  all need to be negated. */
      if( 1 == ((x + y) % 2) )
        theOutput [y] [x] = - theOutput [y] [x];
    }

  return true;
}

И завершите это немного более низким кодом тестирования:

void printMatrix( const double theMatrix [/*Y=*/3] [/*X=*/3] )
{
  for ( int y = 0;  y < 3;  y ++ )
  {
    cout << "[  ";
    for ( int x = 0;  x < 3;  x ++   )
      cout << theMatrix [y] [x] << "  ";
    cout << "]" << endl;
  }
  cout << endl;
}

void matrixMultiply(  const double theMatrixA [/*Y=*/3] [/*X=*/3],
                      const double theMatrixB [/*Y=*/3] [/*X=*/3],
                            double theOutput  [/*Y=*/3] [/*X=*/3]  )
{
  for (   int y = 0;  y < 3;  y ++ )
    for ( int x = 0;  x < 3;  x ++   )
    {
      theOutput [y] [x] = 0;
      for ( int i = 0;  i < 3;  i ++ )
        theOutput [y] [x] +=  theMatrixA [y] [i] * theMatrixB [i] [x];
    }
}

int
main(int argc, char **argv)
{
  if ( argc > 1 )
    SRANDOM( atoi( argv[1] ) );

  double m[3][3] = { { RANDOM_D(0,1e3), RANDOM_D(0,1e3), RANDOM_D(0,1e3) },
                     { RANDOM_D(0,1e3), RANDOM_D(0,1e3), RANDOM_D(0,1e3) },
                     { RANDOM_D(0,1e3), RANDOM_D(0,1e3), RANDOM_D(0,1e3) } };
  double o[3][3], mm[3][3];

  if ( argc <= 2 )
    cout << fixed << setprecision(3);

  printMatrix(m);
  cout << endl << endl;

  SHOW( determinant(m) );
  cout << endl << endl;

  BOUT( inverse(m, o) );
  printMatrix(m);
  printMatrix(o);
  cout << endl << endl;

  matrixMultiply (m, o, mm );
  printMatrix(m);
  printMatrix(o);
  printMatrix(mm);  
  cout << endl << endl;
}

Запоздалая мысль:

9