/* xdet4.h                      freeware                   [[Email Removed]]  */


/* --------------------------------- FUNCTION ------------------------------  */
/* Do    :                                                                    */
/*                                                                            */
/* Call  :                                                                    */
/* Debug :                                                                    */
/* -------------------------------------------------------------------------- */
fraction cofactorF(
pmatrix mA,
int ci,
int cj)
{
     int i;
     int j;
     int n;
     int tci = 0;
     int tcj = 0;
fraction f;

double pMinor [MXR][MXC*TWOCOL];matrix mMinor ={MXR,MXC*TWOCOL,&pMinor [0][0]};

 n = mA->rows;
                mMinor.rows=(n-1);
                mMinor.cols=(n-1)*TWOCOL;

  tci = 0;
  for (i = 0 ; i < mA->rows ; i++)
   {
    if(i != ci){

    tcj = 0;
    for ( j = 0; j < mA->cols ; j++,j++)
     {
      if(j != cj * TWOCOL)
      {
       *(mMinor.pblock+tci *mMinor.cols+tcj   ) = *(mA->pblock+i *mA->cols+j  );
       *(mMinor.pblock+tci *mMinor.cols+ tcj+1) = *(mA->pblock+i *mA->cols+j+1);

       tcj++;
       tcj++;
      }
     }
    tci++; }
   }

          f = detF(&mMinor);
    f.numer = pow(-1,(ci+1+cj+1))* f.numer;

  return(f);
}


/* --------------------------------- FUNCTION ------------------------------  */
/* Do    :                                                                    */
/*                                                                            */
/* Call  :                                                                    */
/* Debug :                                                                    */
/* -------------------------------------------------------------------------- */
void    mcofactorF(
pmatrix mA,
pmatrix mCofactor)
{
     int i;
     int j;
fraction f;

  if (mA->rows != mCofactor->rows || mA-> cols != mCofactor->cols)
    {
     printf("\n mcofactor() error - matrices different sizes");
     printf("\n Press return to continue");
     getchar();
     exit(1);
    }

 for (i = 0 ; i < mA->rows ; i++)
  {
   for ( j = 0 ; j < mA->cols ; j++, j++)
    {
               f =  cofactorF(mA,i, j/TWOCOL);

     *(mCofactor->pblock+i *mCofactor->cols+j  ) = f.numer;
     *(mCofactor->pblock+i *mCofactor->cols+j+1) = f.denom;
    }
  }
}


/* --------------------------------- FUNCTION ------------------------------  */
/* Do    :                                                                    */
/*                                                                            */
/* Call  :                                                                    */
/* Debug :                                                                    */
/* -------------------------------------------------------------------------- */
void    adjointF(
pmatrix mCofactor,
pmatrix mAdjoint)
{
  if (mCofactor->rows != mAdjoint->rows || mCofactor-> cols != mAdjoint->cols)
    {
     printf("\n adjoint() error - matrices different sizes");
     printf("\n Press return to continue");
     getchar();
     exit(1);
    }

 transposemF(mCofactor,mAdjoint);
}


/* --------------------------------- FUNCTION ------------------------------  */
/* Do    :                                                                    */
/*                                                                            */
/* Call  :                                                                    */
/* Debug :                                                                    */
/* -------------------------------------------------------------------------- */
void inverseF(
pmatrix  mA,
pmatrix  mInverse)
{
int      row;
int      n;
fraction f;

double pCofact[MXR][MXC*TWOCOL];matrix mCofact={MXR,MXC*TWOCOL,&pCofact[0][0]};

 n = mA->rows;
               mCofact.rows=n       ;
               mCofact.cols=n*TWOCOL;

     mcofactorF(mA,&mCofact);
       adjointF(   &mCofact,mInverse);

      f = detF(mA);
      f = invF(f);
    for ( row = 0; row < mA->rows ; row++)
     {
        mulrowF(mInverse,row,f);
      frowminiF(mInverse,row  );                        /* a) fraction mini   */
     }
}