Application of the identity matrix II. v05
Submitted By:
xhunga
Rating:
(
Rate It)
/* xbmatop.h freeware [[Email Removed]] */
/* --------------------------------- FUNCTION ------------------------------ */
/* Do : mA + mB = mAplsB */
/* */
/* Call : */
/* Debug : */
/* -------------------------------------------------------------------------- */
void addmF(
pmatrix mA,
pmatrix mB,
pmatrix mAplsB)
{
int i;
int j;
fraction f;
if (mA->rows != mB->rows || mA->cols != mB->cols)
{
printf("\n addm() error - matrices different sizes");
printf("\n Press return to continue");
getchar();
exit(1);
}
if (mA->rows!=mAplsB->rows||mA->cols!=mAplsB->cols)
{
printf("\n addm() 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.numer =
(*(mA->pblock+i *mA->cols+j) *
*(mB->pblock+i *mB->cols+j+1))
+
*(mB->pblock+i *mB->cols+j) *
*(mA->pblock+i *mA->cols+j+1);
f.denom =
*(mA->pblock+i *mA->cols+j+1) *
*(mB->pblock+i *mB->cols+j+1);
f = miniF(f);
*(mAplsB->pblock+i *mAplsB->cols+j) = f.numer;
*(mAplsB->pblock+i *mAplsB->cols+j+1) = f.denom;
}
}
/* --------------------------------- FUNCTION ------------------------------ */
/* Do : mA - mB = mAmnsB */
/* */
/* Call : */
/* Debug : */
/* -------------------------------------------------------------------------- */
void submF(
pmatrix mA,
pmatrix mB,
pmatrix mAsubB)
{
int i;
int j;
fraction f;
if (mA->rows != mB->rows || mA->cols != mB->cols)
{
printf("\n addm() error - matrices different sizes");
printf("\n Press return to continue");
getchar();
exit(1);
}
if (mA->rows!=mAsubB->rows||mA->cols!=mAsubB->cols)
{
printf("\n addm() 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+=2)
{
f.numer =
(*(mA->pblock+i *mA->cols+j) *
*(mB->pblock+i *mB->cols+j+1))
-
*(mB->pblock+i *mB->cols+j) *
*(mA->pblock+i *mA->cols+j+1);
f.denom =
*(mA->pblock+i *mA->cols+j+1) *
*(mB->pblock+i *mB->cols+j+1);
f = miniF(f);
*(mAsubB->pblock+i *mAsubB->cols+j) = f.numer;
*(mAsubB->pblock+i *mAsubB->cols+j+1) = f.denom;
}
}
/* --------------------------------- FUNCTION ------------------------------ */
/* Do : mA * mB = mAB */
/* */
/* Call : */
/* Debug : */
/* -------------------------------------------------------------------------- */
void multmF(
pmatrix mA,
pmatrix mB,
pmatrix mAB)
{
int i;
int j;
int k;
fraction fm;
fraction fa;
for(k= 0; k < mA->rows; k++)
{
for(j= 0; j < mB->cols; j+=2)
{ fa.numer = 0;
fa.denom = 1;
for(i= 0; i < mA->cols; i+=2)
{
fm.numer = *(mA->pblock+k *(mA->cols)+i) *
*(mB->pblock+i/2 *(mB->cols)+j);
fm.denom = *(mA->pblock+k *(mA->cols)+i+1) *
*(mB->pblock+i/2 *(mB->cols)+j+1);
fa = addF(fa,fm);
}
fa = miniF(fa);
*(mAB->pblock+k *(mAB->cols)+j) = fa.numer;
*(mAB->pblock+k *(mAB->cols)+j+1) = fa.denom;
}
}
}
/* --------------------------------- FUNCTION ------------------------------ */
/* Do : s * mA = msA */
/* */
/* Call : */
/* Debug : */
/* -------------------------------------------------------------------------- */
void smultmF(
fraction s,
pmatrix mA,
pmatrix msA)
{
int i;
int j;
fraction f;
for ( i = 0; i < mA->rows ; i++)
for ( j = 0; j < mA->cols ; j+=2)
{
f.numer = (*(mA->pblock+i *mA->cols+j )) * s.numer;
f.denom = (*(mA->pblock+i *mA->cols+j+1)) * s.denom;
f = miniF(f);
(*(msA->pblock+i *msA->cols+j )) = f.numer;
(*(msA->pblock+i *msA->cols+j+1)) = f.denom;
}
}
/* --------------------------------- FUNCTION ------------------------------ */
/* Do : */
/* */
/* Call : */
/* Debug : */
/* -------------------------------------------------------------------------- */
void transposemF(
pmatrix mA,
pmatrix mTrpsA)
{
int i;
int j;
if (mTrpsA->rows * TWOCOL != mA->cols)
{
printf("\n transpose() error - dest matrix rows must = source matrix cols");
printf("\n Press return to continue");
getchar();
exit(1);
}
if (mA->rows * TWOCOL != mTrpsA->cols)
{
printf("\n transpose() error - source matrix rows must = dest matrix cols");
printf("\n Press return to continue");
getchar();
exit(1);
}
for ( i = 0 ; i < mA->rows ; i++)
{
for (j = 0 ; j < mA->cols ; j+=2)
{
*(mTrpsA->pblock+(j/TWOCOL) *mTrpsA->cols+i*TWOCOL) =
*( mA->pblock+i *mA->cols+j);
*(mTrpsA->pblock+(j/TWOCOL) *mTrpsA->cols+i*TWOCOL+1) =
*( mA->pblock+i *mA->cols+j+1);
}
}
}
/* --------------------------------- FUNCTION ------------------------------ */
/* Do : */
/* */
/* Call : */
/* Debug : */
/* -------------------------------------------------------------------------- */
fraction tracemF(
pmatrix m)
{
int i;
fraction f;
fraction fT;
if ((m->rows * TWOCOL) != m->cols)
{
printf("\n trF() error - Square matrix, please");
printf("\n Press return to continue");
getchar();
exit(1);
}
f.numer = 0;
f.denom = 1;
for( i = 0; i < m->rows ; i++)
{
fT.numer = (*(m->pblock+i *m->cols+i*TWOCOL ));
fT.denom = (*(m->pblock+i *m->cols+i*TWOCOL+1));
fT = miniF( fT);
f = addF(f,fT);
f = miniF(f);
}
return(f);
}
/* --------------------------------- FUNCTION ------------------------------ */
/* Do : mA**n = mApn */
/* */
/* Call : */
/* Debug : */
/* -------------------------------------------------------------------------- */
void powermF(
pmatrix mA,
int pn,
pmatrix mApn
)
{
int i;
int n;
double pbT1 [MXR][MXC*TWOCOL];matrix mT1 ={MXR,MXC*TWOCOL,&pbT1 [0][0]};
double pbT2 [MXR][MXC*TWOCOL];matrix mT2 ={MXR,MXC*TWOCOL,&pbT2 [0][0]};
n = mA->rows;
mT1.rows=n; mT1.cols=n*TWOCOL;
mT2.rows=n; mT2.cols=n*TWOCOL;
copymF(mA,mApn);
if(!pn){mid(mApn);}
else
{
copymF(mApn,&mT1);
for(i = pn-1 ; i ; i--)
{
multmF(mApn,&mT1,&mT2);
copymF(&mT2,mApn);
}
}
}
