Вот вариант ответа Бэтти, но он вычисляет правильный обратный. версия Бэтти вычисляет транспонирование обратного.

// computes the inverse of a matrix m
double det = m(0, 0) * (m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2)) -
             m(0, 1) * (m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0)) +
             m(0, 2) * (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0));

double invdet = 1 / det;

Matrix33d minv; // inverse of matrix m
minv(0, 0) = (m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2)) * invdet;
minv(0, 1) = (m(0, 2) * m(2, 1) - m(0, 1) * m(2, 2)) * invdet;
minv(0, 2) = (m(0, 1) * m(1, 2) - m(0, 2) * m(1, 1)) * invdet;
minv(1, 0) = (m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2)) * invdet;
minv(1, 1) = (m(0, 0) * m(2, 2) - m(0, 2) * m(2, 0)) * invdet;
minv(1, 2) = (m(1, 0) * m(0, 2) - m(0, 0) * m(1, 2)) * invdet;
minv(2, 0) = (m(1, 0) * m(2, 1) - m(2, 0) * m(1, 1)) * invdet;
minv(2, 1) = (m(2, 0) * m(0, 1) - m(0, 0) * m(2, 1)) * invdet;
minv(2, 2) = (m(0, 0) * m(1, 1) - m(1, 0) * m(0, 1)) * invdet;

Этот фрагмент кода вычисляет транспонированную обратную матрицу A:

double determinant =    +A(0,0)*(A(1,1)*A(2,2)-A(2,1)*A(1,2))
                        -A(0,1)*(A(1,0)*A(2,2)-A(1,2)*A(2,0))
                        +A(0,2)*(A(1,0)*A(2,1)-A(1,1)*A(2,0));
double invdet = 1/determinant;
result(0,0) =  (A(1,1)*A(2,2)-A(2,1)*A(1,2))*invdet;
result(1,0) = -(A(0,1)*A(2,2)-A(0,2)*A(2,1))*invdet;
result(2,0) =  (A(0,1)*A(1,2)-A(0,2)*A(1,1))*invdet;
result(0,1) = -(A(1,0)*A(2,2)-A(1,2)*A(2,0))*invdet;
result(1,1) =  (A(0,0)*A(2,2)-A(0,2)*A(2,0))*invdet;
result(2,1) = -(A(0,0)*A(1,2)-A(1,0)*A(0,2))*invdet;
result(0,2) =  (A(1,0)*A(2,1)-A(2,0)*A(1,1))*invdet;
result(1,2) = -(A(0,0)*A(2,1)-A(2,0)*A(0,1))*invdet;
result(2,2) =  (A(0,0)*A(1,1)-A(1,0)*A(0,1))*invdet;

Хотя в вопросе оговаривались несингулярные матрицы, вы все равно можете проверить, равен ли определитель нулю (или очень близкому к нулю) и пометить его каким-то образом, чтобы быть в безопасности.

При всем уважении к нашему неизвестному (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;
}

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

Вы также можете захотеть обнаружить очень большие определители, так как ошибки округления повлияют на вашу точность!

Не пытайтесь сделать это самостоятельно, если вы серьезно относитесь к краевым случаям. Таким образом, хотя многие наивные/простые методы теоретически точны, они могут иметь неприятное численное поведение для почти сингулярных матриц. В частности, вы можете получить ошибки отмены/округления, из-за которых вы получите сколь угодно плохие результаты.

«Правильный» способ - это исключение Гаусса с поворотом строки и столбца, так что вы всегда делите на наибольшее оставшееся числовое значение. (Это также стабильно для матриц NxN). Обратите внимание, что поворот строки сам по себе не охватывает все плохие случаи.

Однако IMO, реализующая это правильно и быстро, не стоит вашего времени - используйте хорошо протестированную библиотеку, и есть куча только заголовков.

Я только что создал класс QMatrix. Он использует встроенный вектор > контейнер. QMatrix.h Он использует метод Джордана-Гаусса для вычисления обратной квадратная матрица.

Вы можете использовать его следующим образом:

#include "QMatrix.h"
#include <iostream>

int main(){
QMatrix<double> A(3,3,true);
QMatrix<double> Result = A.inverse()*A; //should give the idendity matrix

std::cout<<A.inverse()<<std::endl;
std::cout<<Result<<std::endl; // for checking
return 0;
}

Обратная функция реализуется следующим образом:

Дан класс со следующими полями:

template<class T> class QMatrix{
public:
int rows, cols;
std::vector<std::vector<T> > A;

обратная() функция:

template<class T> 
QMatrix<T> QMatrix<T>:: inverse(){
Identity<T> Id(rows); //the Identity Matrix as a subclass of QMatrix.
QMatrix<T> Result = *this; // making a copy and transforming it to the Identity matrix
T epsilon = 0.000001;
for(int i=0;i<rows;++i){
    //check if Result(i,i)==0, if true, switch the row with another

    for(int j=i;j<rows;++j){
        if(std::abs(Result(j,j))<epsilon) { //uses Overloading()(int int) to extract element from Result Matrix
            Result.replace_rows(i,j+1); //switches rows i with j+1
        }
        else break;
    }
    // main part, making a triangular matrix
    Id(i)=Id(i)*(1.0/Result(i,i));
    Result(i)=Result(i)*(1.0/Result(i,i));  // Using overloading ()(int) to get a row form the matrix
    for(int j=i+1;j<rows;++j){
        T temp = Result(j,i);
        Result(j) = Result(j) - Result(i)*temp;
        Id(j) = Id(j) - Id(i)*temp; //doing the same operations to the identity matrix
        Result(j,i)=0; //not necessary, but looks nicer than 10^-15
    }
}

// solving a triangular matrix 
for(int i=rows-1;i>0;--i){
    for(int j=i-1;j>=0;--j){
        T temp = Result(j,i);
        Id(j) = Id(j) - temp*Id(i);
        Result(j)=Result(j)-temp*Result(i);
    }
}

return Id;
}

Довольно приятный (как мне кажется) заголовочный файл, содержащий макросы для большинства матричных операций 2x2, 3x3 и 4x4, был доступен с большинством наборов инструментов OpenGL. Не стандартно, но я видел это в разных местах.

Вы можете проверить это здесь. В конце вы найдете обе инверсии 2x2, 3x3 и 4x4.

vvector.h

Я бы также рекомендовал Ilmbase, который является частью OpenEXR. Это хороший набор шаблонных 2,3,4-векторных и матричных подпрограмм.

//Title: Matrix Header File
//Writer: Say OL
//This is a beginner code not an expert one
//No responsibilty for any errors
//Use for your own risk
using namespace std;
int row,col,Row,Col;
double Coefficient;
//Input Matrix
void Input(double Matrix[9][9],int Row,int Col)
{
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
        {
            cout<<"e["<<row<<"]["<<col<<"]=";
            cin>>Matrix[row][col];
        }
}
//Output Matrix
void Output(double Matrix[9][9],int Row,int Col)
{
    for(row=1;row<=Row;row++)
    {
        for(col=1;col<=Col;col++)
            cout<<Matrix[row][col]<<"\t";
        cout<<endl;
    }
}
//Copy Pointer to Matrix
void CopyPointer(double (*Pointer)[9],double Matrix[9][9],int Row,int Col)
{
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
            Matrix[row][col]=Pointer[row][col];
}
//Copy Matrix to Matrix
void CopyMatrix(double MatrixInput[9][9],double MatrixTarget[9][9],int Row,int Col)
{
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
            MatrixTarget[row][col]=MatrixInput[row][col];
}
//Transpose of Matrix
double MatrixTran[9][9];
double (*(Transpose)(double MatrixInput[9][9],int Row,int Col))[9]
{
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
            MatrixTran[col][row]=MatrixInput[row][col];
    return MatrixTran;
}
//Matrix Addition
double MatrixAdd[9][9];
double (*(Addition)(double MatrixA[9][9],double MatrixB[9][9],int Row,int Col))[9]
{
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
            MatrixAdd[row][col]=MatrixA[row][col]+MatrixB[row][col];
    return MatrixAdd;
}
//Matrix Subtraction
double MatrixSub[9][9];
double (*(Subtraction)(double MatrixA[9][9],double MatrixB[9][9],int Row,int Col))[9]
{
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
            MatrixSub[row][col]=MatrixA[row][col]-MatrixB[row][col];
    return MatrixSub;
}
//Matrix Multiplication
int mRow,nCol,pCol,kcol;
double MatrixMult[9][9];
double (*(Multiplication)(double MatrixA[9][9],double MatrixB[9][9],int mRow,int nCol,int pCol))[9]
{
    for(row=1;row<=mRow;row++)
        for(col=1;col<=pCol;col++)
        {
            MatrixMult[row][col]=0.0;
            for(kcol=1;kcol<=nCol;kcol++)
                MatrixMult[row][col]+=MatrixA[row][kcol]*MatrixB[kcol][col];
        }
    return MatrixMult;
}
//Interchange Two Rows
double RowTemp[9][9];
double MatrixInter[9][9];
double (*(InterchangeRow)(double MatrixInput[9][9],int Row,int Col,int iRow,int jRow))[9]
{
    CopyMatrix(MatrixInput,MatrixInter,Row,Col);
    for(col=1;col<=Col;col++)
    {
        RowTemp[iRow][col]=MatrixInter[iRow][col];
        MatrixInter[iRow][col]=MatrixInter[jRow][col];
        MatrixInter[jRow][col]=RowTemp[iRow][col];
    }
    return MatrixInter;
}
//Pivote Downward
double MatrixDown[9][9];
double (*(PivoteDown)(double MatrixInput[9][9],int Row,int Col,int tRow,int tCol))[9]
{
    CopyMatrix(MatrixInput,MatrixDown,Row,Col);
    Coefficient=MatrixDown[tRow][tCol];
    if(Coefficient!=1.0)
        for(col=1;col<=Col;col++)
            MatrixDown[tRow][col]/=Coefficient;
    if(tRow<Row)
        for(row=tRow+1;row<=Row;row++)
        {
            Coefficient=MatrixDown[row][tCol];
            for(col=1;col<=Col;col++)
                MatrixDown[row][col]-=Coefficient*MatrixDown[tRow][col];
        }
return MatrixDown;
}
//Pivote Upward
double MatrixUp[9][9];
double (*(PivoteUp)(double MatrixInput[9][9],int Row,int Col,int tRow,int tCol))[9]
{
    CopyMatrix(MatrixInput,MatrixUp,Row,Col);
    Coefficient=MatrixUp[tRow][tCol];
    if(Coefficient!=1.0)
        for(col=1;col<=Col;col++)
            MatrixUp[tRow][col]/=Coefficient;
    if(tRow>1)
        for(row=tRow-1;row>=1;row--)
        {
            Coefficient=MatrixUp[row][tCol];
            for(col=1;col<=Col;col++)
                MatrixUp[row][col]-=Coefficient*MatrixUp[tRow][col];
        }
    return MatrixUp;
}
//Pivote in Determinant
double MatrixPiv[9][9];
double (*(Pivote)(double MatrixInput[9][9],int Dim,int pTarget))[9]
{
    CopyMatrix(MatrixInput,MatrixPiv,Dim,Dim);
    for(row=pTarget+1;row<=Dim;row++)
    {
        Coefficient=MatrixPiv[row][pTarget]/MatrixPiv[pTarget][pTarget];
        for(col=1;col<=Dim;col++)
        {
            MatrixPiv[row][col]-=Coefficient*MatrixPiv[pTarget][col];
        }
    }
    return MatrixPiv;
}
//Determinant of Square Matrix
int dCounter,dRow;
double Det;
double MatrixDet[9][9];
double Determinant(double MatrixInput[9][9],int Dim)
{
    CopyMatrix(MatrixInput,MatrixDet,Dim,Dim);
    Det=1.0;
    if(Dim>1)
    {
        for(dRow=1;dRow<Dim;dRow++)
        {
            dCounter=dRow;
            while((MatrixDet[dRow][dRow]==0.0)&(dCounter<=Dim))
            {
                dCounter++;
                Det*=-1.0;
                CopyPointer(InterchangeRow(MatrixDet,Dim,Dim,dRow,dCounter),MatrixDet,Dim,Dim);
            }
            if(MatrixDet[dRow][dRow]==0)
            {
                Det=0.0;
                break;
            }
            else
            {
                Det*=MatrixDet[dRow][dRow];
                CopyPointer(Pivote(MatrixDet,Dim,dRow),MatrixDet,Dim,Dim);
            }
        }
        Det*=MatrixDet[Dim][Dim];
    }
    else Det=MatrixDet[1][1];
    return Det;
}
//Matrix Identity
double MatrixIdent[9][9];
double (*(Identity)(int Dim))[9]
{
    for(row=1;row<=Dim;row++)
        for(col=1;col<=Dim;col++)
            if(row==col)
                MatrixIdent[row][col]=1.0;
            else
                MatrixIdent[row][col]=0.0;
    return MatrixIdent;
}
//Join Matrix to be Augmented Matrix
double MatrixJoin[9][9];
double (*(JoinMatrix)(double MatrixA[9][9],double MatrixB[9][9],int Row,int ColA,int ColB))[9]
{
    Col=ColA+ColB;
    for(row=1;row<=Row;row++)
        for(col=1;col<=Col;col++)
            if(col<=ColA)
                MatrixJoin[row][col]=MatrixA[row][col];
            else
                MatrixJoin[row][col]=MatrixB[row][col-ColA];
    return MatrixJoin;
}
//Inverse of Matrix
double (*Pointer)[9];
double IdentMatrix[9][9];
int Counter;
double MatrixAug[9][9];
double MatrixInv[9][9];
double (*(Inverse)(double MatrixInput[9][9],int Dim))[9]
{
    Row=Dim;
    Col=Dim+Dim;
    Pointer=Identity(Dim);
    CopyPointer(Pointer,IdentMatrix,Dim,Dim);
    Pointer=JoinMatrix(MatrixInput,IdentMatrix,Dim,Dim,Dim);
    CopyPointer(Pointer,MatrixAug,Row,Col);
    for(Counter=1;Counter<=Dim;Counter++)   
    {
        Pointer=PivoteDown(MatrixAug,Row,Col,Counter,Counter);
        CopyPointer(Pointer,MatrixAug,Row,Col);
    }
    for(Counter=Dim;Counter>1;Counter--)
    {
        Pointer=PivoteUp(MatrixAug,Row,Col,Counter,Counter);
        CopyPointer(Pointer,MatrixAug,Row,Col);
    }
    for(row=1;row<=Dim;row++)
        for(col=1;col<=Dim;col++)
            MatrixInv[row][col]=MatrixAug[row][col+Dim];
    return MatrixInv;
}
//Gauss-Jordan Elemination
double MatrixGJ[9][9];
double VectorGJ[9][9];
double (*(GaussJordan)(double MatrixInput[9][9],double VectorInput[9][9],int Dim))[9]
{
    Row=Dim;
    Col=Dim+1;
    Pointer=JoinMatrix(MatrixInput,VectorInput,Dim,Dim,1);
    CopyPointer(Pointer,MatrixGJ,Row,Col);
    for(Counter=1;Counter<=Dim;Counter++)   
    {
        Pointer=PivoteDown(MatrixGJ,Row,Col,Counter,Counter);
        CopyPointer(Pointer,MatrixGJ,Row,Col);
    }
    for(Counter=Dim;Counter>1;Counter--)
    {
        Pointer=PivoteUp(MatrixGJ,Row,Col,Counter,Counter);
        CopyPointer(Pointer,MatrixGJ,Row,Col);
    }
    for(row=1;row<=Dim;row++)
        for(col=1;col<=1;col++)
            VectorGJ[row][col]=MatrixGJ[row][col+Dim];
    return VectorGJ;
}
//Generalized Gauss-Jordan Elemination
double MatrixGGJ[9][9];
double VectorGGJ[9][9];
double (*(GeneralizedGaussJordan)(double MatrixInput[9][9],double VectorInput[9][9],int Dim,int vCol))[9]
{
    Row=Dim;
    Col=Dim+vCol;
    Pointer=JoinMatrix(MatrixInput,VectorInput,Dim,Dim,vCol);
    CopyPointer(Pointer,MatrixGGJ,Row,Col);
    for(Counter=1;Counter<=Dim;Counter++)   
    {
        Pointer=PivoteDown(MatrixGGJ,Row,Col,Counter,Counter);
        CopyPointer(Pointer,MatrixGGJ,Row,Col);
    }
    for(Counter=Dim;Counter>1;Counter--)
    {
        Pointer=PivoteUp(MatrixGGJ,Row,Col,Counter,Counter);
        CopyPointer(Pointer,MatrixGGJ,Row,Col);
    }
    for(row=1;row<=Row;row++)
        for(col=1;col<=vCol;col++)
            VectorGGJ[row][col]=MatrixGGJ[row][col+Dim];
    return VectorGGJ;
}
//Matrix Sparse, Three Diagonal Non-Zero Elements
double MatrixSpa[9][9];
double (*(Sparse)(int Dimension,double FirstElement,double SecondElement,double ThirdElement))[9]
{
    MatrixSpa[1][1]=SecondElement;
    MatrixSpa[1][2]=ThirdElement;
    MatrixSpa[Dimension][Dimension-1]=FirstElement;
    MatrixSpa[Dimension][Dimension]=SecondElement;
    for(int Counter=2;Counter<Dimension;Counter++)
    {
        MatrixSpa[Counter][Counter-1]=FirstElement;
        MatrixSpa[Counter][Counter]=SecondElement;
        MatrixSpa[Counter][Counter+1]=ThirdElement;
    }
    return MatrixSpa;
}

Скопируйте и сохраните приведенный выше код как Matrix.h, затем попробуйте следующий код:

#include<iostream>
#include<conio.h>
#include"Matrix.h"
int Dim;
double Matrix[9][9];
int main()
{
    cout<<"Enter your matrix dimension: ";
    cin>>Dim;
    Input(Matrix,Dim,Dim);
    cout<<"Your matrix:"<<endl;
    Output(Matrix,Dim,Dim);
    cout<<"The inverse:"<<endl;
    Output(Inverse(Matrix,Dim),Dim,Dim);
    getch();
}

Я пошел дальше и написал его на питоне, так как я думаю, что он намного читабельнее, чем на С++ для такой проблемы. Порядок функций соответствует порядку операций для решения этого вручную через это видео. Просто импортируйте это и вызовите "print_invert" на своей матрице.

def print_invert (matrix):
  i_matrix = invert_matrix (matrix)
  for line in i_matrix:
    print (line)
  return

def invert_matrix (matrix):
  determinant = str (determinant_of_3x3 (matrix))
  cofactor = make_cofactor (matrix)
  trans_matrix = transpose_cofactor (cofactor)

  trans_matrix[:] = [[str (element) +'/'+ determinant for element in row] for row in trans_matrix]

  return trans_matrix

def determinant_of_3x3 (matrix):
  multiplication = 1
  neg_multiplication = 1
  total = 0
  for start_column in range (3):
    for row in range (3):
      multiplication *= matrix[row][(start_column+row)%3]
      neg_multiplication *= matrix[row][(start_column-row)%3]
    total += multiplication - neg_multiplication
    multiplication = neg_multiplication = 1
  if total == 0:
    total = 1
  return total

def make_cofactor (matrix):
  cofactor = [[0,0,0],[0,0,0],[0,0,0]]
  matrix_2x2 = [[0,0],[0,0]]
  # For each element in matrix...
  for row in range (3):
    for column in range (3):

      # ...make the 2x2 matrix in this inner loop
      matrix_2x2 = make_2x2_from_spot_in_3x3 (row, column, matrix)
      cofactor[row][column] = determinant_of_2x2 (matrix_2x2)

  return flip_signs (cofactor)

def make_2x2_from_spot_in_3x3 (row, column, matrix):
  c_count = 0
  r_count = 0
  matrix_2x2 = [[0,0],[0,0]]
  # ...make the 2x2 matrix in this inner loop
  for inner_row in range (3):
    for inner_column in range (3):
      if row is not inner_row and inner_column is not column:
        matrix_2x2[r_count % 2][c_count % 2] = matrix[inner_row][inner_column]
        c_count += 1
    if row is not inner_row:
      r_count += 1
  return matrix_2x2

def determinant_of_2x2 (matrix):
  total = matrix[0][0] * matrix [1][1]
  return total - (matrix [1][0] * matrix [0][1])

def flip_signs (cofactor):
  sign_pos = True 
  # For each element in matrix...
  for row in range (3):
    for column in range (3):
      if sign_pos:
        sign_pos = False
      else:
        cofactor[row][column] *= -1
        sign_pos = True
  return cofactor

def transpose_cofactor (cofactor):
  new_cofactor = [[0,0,0],[0,0,0],[0,0,0]]
  for row in range (3):
    for column in range (3):
      new_cofactor[column][row] = cofactor[row][column]
  return new_cofactor