matrix.c

Go to the documentation of this file.

//DEBUG

/*
 * \file
 * This file is in the public domain.
 */

#include "matrix.h"

typedef enum {
 MatrixIdentity=0x01,
 MatrixTranslation=0x02,
 MatrixScaling=0x04,
 MatrixUnknown=0x08
} MATRIXFLAGS;

typedef struct{
 double eM11;
 double eM12;
 double eM21;
 double eM22;
 double eDx;
 double eDy;
 DWORD flags;
 DWORD padding;
} MATRIXINTERNAL;


#ifdef DEBUG
#define ASSERTMATRIX(x) AssertMatrix(x,__FILE__,__LINE__)
#else
#define ASSERTMATRIX(x)
#endif
static void AssertMatrix(MATRIX *matrix, char *file, int line){
 DWORD flags=0;
 MATRIXINTERNAL *mat=(MATRIXINTERNAL*)matrix;
 if(!mat)return;
 if(mat->eM12!=0.0 || mat->eM21!=0.0){
  flags|=MatrixUnknown;
 } else if(mat->eM11==1.0 && mat->eM22==1.0){
  flags|=MatrixIdentity;
 } else {
  flags|=MatrixScaling;
 }
 if(mat->eDx!=0.0 || mat->eDy!=0.0){
  flags|=MatrixTranslation;
 }
 if((mat->flags&MatrixIdentity)!=(flags&MatrixIdentity)||
    ((flags&MatrixScaling)&&!(mat->flags&(MatrixScaling|MatrixIdentity)))||
    ((flags&MatrixTranslation)&&!(mat->flags&MatrixTranslation))
  ){
//  DebugOut("Matrix assertion failed: calcflags: %d; matrix flags: %d [%s line %d]",flags,mat->flags,file,line);
//  DebugOut("Matrix: %f %f %f %f %f %f",mat->eM11,mat->eM12,mat->eM21,mat->eM22,mat->eDx,mat->eDy);
//  ASSERT(FALSE);
 }
}

BOOL MatrixEquals(MATRIX *matrix, MATRIX *matrix2){
 MATRIXINTERNAL *mat1=(MATRIXINTERNAL*)matrix;
 MATRIXINTERNAL *mat2=(MATRIXINTERNAL*)matrix2;
 if(mat1==mat2)return TRUE;
 if(!mat1||!mat2)return FALSE;
 return (mat1->eM11==mat2->eM11&&
  mat1->eM12==mat2->eM12&&
  mat1->eM21==mat2->eM21&&
  mat1->eM22==mat2->eM22&&
  mat1->eDx==mat2->eDx&&
  mat1->eDy==mat2->eDy);
}


BOOL MatrixAreClose(MATRIX *matrix, MATRIX *matrix2){
 MATRIXINTERNAL *mat1=(MATRIXINTERNAL*)matrix;
 MATRIXINTERNAL *mat2=(MATRIXINTERNAL*)matrix2;
 if(mat1==mat2)return TRUE;
 if(!mat1||!mat2)return FALSE;
 return (abs(mat1->eM11-mat2->eM11)<=1e-12&&
  abs(mat1->eM12-mat2->eM12)<=1e-12&&
  abs(mat1->eM21-mat2->eM21)<=1e-12&&
  abs(mat1->eM22-mat2->eM22)<=1e-12&&
  abs(mat1->eDx-mat2->eDx)<=1e-12&&
  abs(mat1->eDy-mat2->eDy)<=1e-12);
}


void MatrixTransformPoints(MATRIX *matrix, double *pts, DWORD numpts){
 DWORD i;
 MATRIXINTERNAL *mat=(MATRIXINTERNAL*)matrix;
 if(!mat||mat->flags==MatrixIdentity)return;
 for(i=0;i<numpts;i++){
  if(mat->flags==MatrixTranslation){
   pts[0]+=mat->eDx;
   pts[1]+=mat->eDy;
  } else if(mat->flags==MatrixScaling){
   pts[0]*=mat->eM11;
   pts[1]*=mat->eM22;
  } else if(mat->flags==(MatrixScaling|MatrixTranslation)){
   pts[0]*=mat->eM11;
   pts[1]*=mat->eM22;
   pts[0]+=mat->eDx;
   pts[1]+=mat->eDy;
  } else {
   double x,y;
   x=mat->eM11*pts[0]+mat->eM21*pts[1]+mat->eDx;
   y=mat->eM12*pts[0]+mat->eM22*pts[1]+mat->eDy;
   pts[0]=x;
   pts[1]=y;
  }
  pts+=2;
 }
}

void MatrixTransformVectors(MATRIX *matrix, double *pts, DWORD numpts){
 DWORD i;
 MATRIXINTERNAL *mat=(MATRIXINTERNAL*)matrix;
 if(!mat||(mat->flags&MatrixIdentity))return;
 for(i=0;i<numpts;i++){
  if(mat->flags&MatrixScaling){
   pts[0]*=mat->eM11;
   pts[1]*=mat->eM22;
  } else {
   double x,y;
   x=mat->eM11*pts[0]+mat->eM21*pts[1];
   y=mat->eM12*pts[0]+mat->eM22*pts[1];
   pts[0]=x;
   pts[1]=y;
  }
  pts+=2;
 }
}


double MatrixDeterminant(MATRIX *matrix){
 MATRIXINTERNAL *mat=(MATRIXINTERNAL*)matrix;
 if(!mat)return 0.0;
 ASSERTMATRIX(mat);
 if(mat->flags&MatrixIdentity){
  return 1.0;
 } else if(mat->flags&MatrixScaling){
  return mat->eM11*mat->eM22;
 } else {
  return mat->eM11*mat->eM22-mat->eM12*mat->eM21;
 }
}

static inline void MatrixUpdate(MATRIX *matrix){
 MATRIXINTERNAL *mat=(MATRIXINTERNAL*)matrix;
 mat->flags=0;
 if(mat->eM12!=0.0 || mat->eM21!=0.0){
  mat->flags|=MatrixUnknown;
 } else if(mat->eM11==1.0 && mat->eM22==1.0){
  mat->flags|=MatrixIdentity;
 } else {
  mat->flags|=MatrixScaling;
 }
 if(mat->eDx!=0.0 || mat->eDy!=0.0){
  mat->flags|=MatrixTranslation;
 }
}

double MatrixGetValue(MATRIX *matrix, MATRIXELEMENT elem){
 double *values=(double*)matrix;
 if(!values||elem<0||elem>=6)return 0.0;
 return values[elem];
}

BOOL MatrixSetValue(MATRIX *matrix, MATRIXELEMENT elem, double value){
 double *values=(double*)matrix;
 if(!values||elem<0||elem>=6)return FALSE;
 values[elem]=value;
 MatrixUpdate(matrix);
 return TRUE;
}

BOOL MatrixGetValueArray(MATRIX *matrix, double *values){
 MATRIXINTERNAL *mat=(MATRIXINTERNAL *)matrix;
 if(!matrix||!values)return FALSE;
 values[0]=mat->eM11;
 values[1]=mat->eM12;
 values[2]=mat->eM21;
 values[3]=mat->eM22;
 values[4]=mat->eDx;
 values[5]=mat->eDy;
 return TRUE;
}

BOOL MatrixSetValueArray(MATRIX *matrix, double *values){
 MATRIXINTERNAL *mat=(MATRIXINTERNAL *)matrix;
 if(!matrix||!values)return FALSE;
 mat->eM11=values[0];
 mat->eM12=values[1];
 mat->eM21=values[2];
 mat->eM22=values[3];
 mat->eDx=values[4];
 mat->eDy=values[5];
 MatrixUpdate(mat);
 return TRUE;
}

BOOL MatrixIsInvertible(MATRIX *matrix){
 double det;
 det=MatrixDeterminant(matrix);
 return !(det>-1e-12&&det<1e-12);
}

BOOL MatrixInvert(MATRIX *matrixDst, MATRIX *matrixSrc){
 double det;
 MATRIXINTERNAL tmp;
 MATRIXINTERNAL *matDst=(MATRIXINTERNAL*)matrixDst;
 MATRIXINTERNAL *matSrc=(MATRIXINTERNAL*)matrixSrc;
 ASSERTMATRIX(matDst);
 ASSERTMATRIX(matSrc);
 if(matSrc->flags==MatrixIdentity){
  MatrixSetIdentity(matDst);
  return TRUE;
 }
 det=MatrixDeterminant(matSrc);
 if(det>-1e-12&&det<1e-12){
  MatrixSetIdentity(matDst);
  return FALSE;
 }
 tmp.eM11=matSrc->eM22/det;
 tmp.eM12=-matSrc->eM12/det;
 tmp.eM21=-matSrc->eM21/det;
 tmp.eM22=matSrc->eM11/det;
 tmp.eDx=-matSrc->eDx*tmp.eM11-matSrc->eDy*tmp.eM21;
 tmp.eDy=-matSrc->eDx*tmp.eM12-matSrc->eDy*tmp.eM22;
 MatrixUpdate(&tmp);
 *(MATRIXINTERNAL *)matDst=tmp;
 return TRUE;
}

BOOL MatrixIsIdentity(MATRIX *matrix){
 MATRIXINTERNAL *mat=(MATRIXINTERNAL*)matrix;
 if(!mat)return FALSE;
 ASSERTMATRIX(mat);
 return (mat->flags==MatrixIdentity||(
  mat->eM11==1.0&&mat->eM22==1.0&&mat->eM12==0.0&&
  mat->eM21==0.0&&mat->eDx==0.0&&mat->eDy==0.0
 ));
}

void MatrixSetIdentity(MATRIX *matrix){
 MATRIXINTERNAL *mat=(MATRIXINTERNAL*)matrix;
 mat->eM11=1.0;
 mat->eM12=0.0;
 mat->eM21=0.0;
 mat->eM22=1.0;
 mat->eDx=0.0;
 mat->eDy=0.0;
 mat->flags=MatrixIdentity;
}


void MatrixSetValues(MATRIX *matrix, double m11, double m12, double m21, double m22, double dx, double dy){
 MATRIXINTERNAL *mat=(MATRIXINTERNAL*)matrix;
 mat->eM11=m11;
 mat->eM12=m12;
 mat->eM21=m21;
 mat->eM22=m22;
 mat->eDx=dx;
 mat->eDy=dy;
 MatrixUpdate(mat);
}

static void DebugOutMatrixInternal(MATRIXINTERNAL *mat){
 DebugOut("%f %f %f %f %f %f: flags=%04X",mat->eM11,mat->eM12,mat->eM21,mat->eM22,mat->eDx,mat->eDy,mat->flags);
}

void MatrixMultiply(MATRIX *matrixResult, MATRIX *matrixLeft, MATRIX *matrixRight){
 MATRIXINTERNAL mat;
 MATRIXINTERNAL *matResult=(MATRIXINTERNAL*)matrixResult;
 MATRIXINTERNAL *matRight=(MATRIXINTERNAL*)matrixRight;
 MATRIXINTERNAL *matLeft=(MATRIXINTERNAL*)matrixLeft;
 if(!matResult||!matRight||!matLeft)
  return;
// matRight->flags=MatrixUnknown|MatrixTranslation;
// matLeft->flags=MatrixUnknown|MatrixTranslation;
 ASSERTMATRIX(matRight);
 ASSERTMATRIX(matLeft);
 if(matLeft->flags!=MatrixIdentity){
  if(matRight->flags==MatrixIdentity){
    mat=*matLeft;
  } else if(matRight->flags&MatrixIdentity){
    mat.eM11=matLeft->eM11;
    mat.eM12=matLeft->eM12;
    mat.eM21=matLeft->eM21;
    mat.eM22=matLeft->eM22;
    mat.eDx=matLeft->eDx+matRight->eDx;
    mat.eDy=matLeft->eDy+matRight->eDy;
    MatrixUpdate(&mat);
  } else if(matLeft->flags&MatrixIdentity){
    mat.eM11=matRight->eM11;
    mat.eM12=matRight->eM12;
    mat.eM21=matRight->eM21;
    mat.eM22=matRight->eM22;
    mat.eDx=matLeft->eDx*matRight->eM11+matLeft->eDy*matRight->eM21+matRight->eDx;
    mat.eDy=matLeft->eDx*matRight->eM12+matLeft->eDy*matRight->eM22+matRight->eDy;
    MatrixUpdate(&mat);
  } else if(matRight->flags==MatrixScaling && matLeft->flags==MatrixScaling){
    mat.eM21=mat.eM12=mat.eDx=mat.eDy=0.0;
    mat.eM11=matRight->eM11*matLeft->eM11;
    mat.eM22=matRight->eM22*matLeft->eM22;
    mat.flags=MatrixScaling;
  } else if(matRight->flags&MatrixScaling){
   mat.eM11=matRight->eM11*matLeft->eM11;
   mat.eM21=matRight->eM11*matLeft->eM21;
   mat.eM12=matRight->eM22*matLeft->eM12;
   mat.eM22=matRight->eM22*matLeft->eM22;
   if(matLeft->flags&MatrixTranslation){
    mat.eDx=matLeft->eDx*matRight->eM11+matRight->eDx;
    mat.eDy=matLeft->eDy*matRight->eM22+matRight->eDy;
   } else {
    mat.eDx=matRight->eDx;
    mat.eDy=matRight->eDy;
   }
   MatrixUpdate(&mat);
  } else if(matLeft->flags&MatrixScaling){
   mat.eM11=matRight->eM11*matLeft->eM11;
   mat.eM21=matRight->eM21*matLeft->eM22;
   mat.eM12=matRight->eM12*matLeft->eM11;
   mat.eM22=matRight->eM22*matLeft->eM22;
   if(matLeft->flags&MatrixTranslation){
    mat.eDx=matLeft->eDx*matRight->eM11+matLeft->eDy*matRight->eM21+matRight->eDx;
    mat.eDy=matLeft->eDx*matRight->eM12+matLeft->eDy*matRight->eM22+matRight->eDy;
   } else {
    mat.eDx=matRight->eDx;
    mat.eDy=matRight->eDy;
   }
   MatrixUpdate(&mat);
  } else {
   mat.eM11=matRight->eM11*matLeft->eM11+matRight->eM21*matLeft->eM12;
   mat.eM21=matRight->eM11*matLeft->eM21+matRight->eM21*matLeft->eM22;
   mat.eM12=matRight->eM12*matLeft->eM11+matRight->eM22*matLeft->eM12;
   mat.eM22=matRight->eM12*matLeft->eM21+matRight->eM22*matLeft->eM22;
   if(matLeft->flags&MatrixTranslation){
    mat.eDx=matLeft->eDx*matRight->eM11+matLeft->eDy*matRight->eM21+matRight->eDx;
    mat.eDy=matLeft->eDx*matRight->eM12+matLeft->eDy*matRight->eM22+matRight->eDy;
   } else {
    mat.eDx=matRight->eDx;
    mat.eDy=matRight->eDy;
   }
   MatrixUpdate(&mat);
  }
 } else {
  mat=*matRight;
 }
 *matResult=mat;
}

static void MatrixApply(MATRIX *matrix, MATRIXINTERNAL *matSrc, MATRIXMODE mode){
 if(mode==MatrixSet || ((MATRIXINTERNAL*)matrix)->flags==MatrixIdentity){
  *(MATRIXINTERNAL*)matrix=*matSrc;
 } else if(mode==MatrixAppend){
  MatrixMultiply(matrix,matrix,(MATRIX*)matSrc);
 } else if(mode==MatrixPrepend){
  MatrixMultiply(matrix,(MATRIX*)matSrc,matrix);
 }
}

void MatrixTranslate(MATRIX *matrix, double x, double y, MATRIXMODE mode){
 MATRIXINTERNAL *mat=(MATRIXINTERNAL*)matrix;
 ASSERTMATRIX(mat);
 if(mode==MatrixSet || mat->flags==MatrixIdentity){
  mat->eM11=1.0;
  mat->eM12=0.0;
  mat->eM21=0.0;
  mat->eM22=1.0;
  mat->eDx=x;
  mat->eDy=y;
  mat->flags=MatrixTranslation|MatrixIdentity;
 } else if(mode==MatrixAppend){
  mat->eDx+=x;
  mat->eDy+=y;
  mat->flags|=MatrixTranslation;
 } else if(mode==MatrixPrepend){
  if(mat->flags&MatrixIdentity){
   mat->eDx+=x;
   mat->eDy+=y;
  } else if(mat->flags&MatrixScaling){
   mat->eDx+=x*mat->eM11;
   mat->eDy+=y*mat->eM22;
  } else {
   mat->eDx+=x*mat->eM11+y*mat->eM21;
   mat->eDy+=x*mat->eM12+y*mat->eM22;
  }
  mat->flags|=MatrixTranslation;
 }
}

void MatrixRotate(MATRIX *matrix, double degrees, MATRIXMODE mode){
 MATRIXINTERNAL mat;
 double rad,sinrad,cosrad;
 if(!matrix)return;
 ASSERTMATRIX(matrix);
 rad=degrees*DEGTORAD;
 sinrad=sin(rad);
 cosrad=cos(rad);
 mat.eM11=cosrad;
 mat.eM12=sinrad;
 mat.eM21=-sinrad;
 mat.eM22=cosrad;
 mat.eDx=0.0;
 mat.eDy=0.0;
 mat.flags=MatrixUnknown;
 MatrixApply(matrix,&mat,mode);
}

void MatrixRotateAt(MATRIX *matrix, double degrees, double x, double y, MATRIXMODE mode){
 MATRIXINTERNAL mat;
 double rad,sinrad,cosrad;
 ASSERTMATRIX(matrix);
 rad=degrees*DEGTORAD;
 sinrad=sin(rad);
 cosrad=cos(rad);
 mat.eM11=cosrad;
 mat.eM12=sinrad;
 mat.eM21=-sinrad;
 mat.eM22=cosrad;
 mat.eDx=-x*cosrad+y*sinrad+x;
 mat.eDy=-x*sinrad-y*cosrad+y;
 mat.flags=MatrixUnknown|MatrixTranslation;
 MatrixApply(matrix,&mat,mode);
}


void MatrixScale(MATRIX *matrix, double x, double y, MATRIXMODE mode){
 MATRIXINTERNAL mat;
 ASSERTMATRIX(matrix);
 mat.eM11=x;
 mat.eM12=0.0;
 mat.eM21=0.0;
 mat.eM22=y;
 mat.eDx=0.0;
 mat.eDy=0.0;
 mat.flags=MatrixScaling;
 MatrixApply(matrix,&mat,mode);
}

void MatrixScaleAt(MATRIX *matrix, double x, double y, double centerX, double centerY, MATRIXMODE mode){
 MATRIXINTERNAL mat;
 ASSERTMATRIX(matrix);
 mat.eM11=x;
 mat.eM12=0.0;
 mat.eM21=0.0;
 mat.eM22=y;
 mat.eDx=centerX-(x*centerX);
 mat.eDy=centerY-(y*centerY);
 mat.flags=MatrixScaling|MatrixTranslation;
 MatrixApply(matrix,&mat,mode);
}

void MatrixSkew(MATRIX *matrix, double angleX, double angleY, MATRIXMODE mode){
 MATRIXINTERNAL mat;
 if(!matrix)return;
 ASSERTMATRIX(matrix);
 mat.eM11=1.0;
 mat.eM12=tan(angleY*DEGTORAD);
 mat.eM21=tan(angleX*DEGTORAD);
 mat.eM22=1.0;
 mat.eDx=0.0;
 mat.eDy=0.0;
 mat.flags=MatrixUnknown;
 MatrixApply(matrix,&mat,mode);
}


void MatrixShear(MATRIX *matrix, double shearX, double shearY, MATRIXMODE mode){
 MATRIXINTERNAL mat;
 if(!matrix)return;
 ASSERTMATRIX(matrix);
 mat.eM11=1.0;
 mat.eM12=shearY;
 mat.eM21=shearX;
 mat.eM22=1.0;
 mat.eDx=0.0;
 mat.eDy=0.0;
 mat.flags=MatrixUnknown;
 MatrixApply(matrix,&mat,mode);
}

void MatrixParallelogram(MATRIX *matrix, double *rect, double *points, MATRIXMODE mode){
 MATRIXINTERNAL mat;
 if(!matrix||!rect||!points)return;
 ASSERTMATRIX(matrix);
 if(rect[2]==0 || rect[3]==0){
  mat.eM11=mat.eM22=1.0;
  mat.eM12=mat.eM21=0.0;  
  mat.eDx=mat.eDx=0.0;
  mat.flags=MatrixIdentity;
 } else {
  mat.eM11 = (points[2] - points[0]) / rect[2];
  mat.eM12 = (points[3] - points[1]) / rect[2];
  mat.eM21 = (points[4] - points[0]) / rect[3];
  mat.eM22 = (points[5] - points[1]) / rect[3];
  mat.eDx=-rect[0]*mat.eM11-rect[1]*mat.eM21+points[0];
  mat.eDy=-rect[0]*mat.eM12-rect[1]*mat.eM22+points[1];
  mat.flags=MatrixUnknown|MatrixTranslation;
 }
 MatrixApply(matrix,&mat,mode);
}

#if 0
#define ASSERT2(x,file,line) if(!(x)){\
  DebugOut("Assertion failed (%s line %d): %s",file,line,#x); }

#define AME(m,a,b,c,d,e,f) \
  AssertMatEqual(m,a,b,c,d,e,f,__FILE__,__LINE__)

static BOOL AreCloseAdv(double value1, double value2){
 if (value1 == value2)return TRUE;
 double num = ((abs(value1) + abs(value2)) + 10.0) * 2.2204460492503131E-16;
 double num2 = value1 - value2;
 return ((-num < num2) && (num > num2));
}
static BOOL AreClose(double x, double y){
 return abs(x-y)<=0.0001;
}

static void AssertMatEqual(MATRIX *matrix, 
  double m11, double m12, double m21, double m22, double dx, double dy,
  char *file, int line){
 ASSERT2(AreClose(GetM11(matrix),m11),file,line); 
 ASSERT2(AreClose(GetM12(matrix),m12),file,line); 
 ASSERT2(AreClose(GetM21(matrix),m21),file,line); 
 ASSERT2(AreClose(GetM22(matrix),m22),file,line); 
 ASSERT2(AreClose(GetDx(matrix),dx),file,line); 
 ASSERT2(AreClose(GetDy(matrix),dy),file,line); 
}
/*
int main(int argc,char**argv){
 MATRIX m,m2;
 DebugOut("Start of tests");
 MatrixSetIdentity(&m);
 AME(&m,1,0,0,1,0,0);
 ASSERT(MatrixIsIdentity(&m));
 MatrixSetValues(&m,123, 24, 82, 16, 47, 30);
 ASSERT(!MatrixIsInvertible(&m));
 ASSERT(!MatrixIsIdentity(&m));
 ASSERT(GetDx(&m)==47);
 ASSERT(GetDy(&m)==30);
 MatrixSetValues(&m,156, 46, 0, 0, 106, 19);
 ASSERT(!MatrixIsInvertible(&m));
 MatrixSetValues(&m,146, 66, 158, 104, 42, 150);
 ASSERT(MatrixIsInvertible(&m));
 MatrixSetValues(&m,119, 140, 145, 74, 102, 58);
 ASSERT(MatrixIsInvertible(&m));
 MatrixSetIdentity(&m);
 AME(&m,1,0,0,1,0,0);
 MatrixSetValues(&m,10,20,30,40,50,60);
 MatrixRotate(&m,180,MatrixPrepend);
 AME(&m,-10,-20,-30,-40,50,60);
 MatrixSetIdentity(&m);
 ASSERT(MatrixIsIdentity(&m));
 MatrixRotate(&m,45,MatrixPrepend);
 ASSERT(!MatrixIsIdentity(&m));
 AME(&m,0.7071068,0.7071068,-0.7071068,0.7071068,0,0);
 MatrixRotate(&m,135,MatrixPrepend);
 AME(&m,-1,0,0,-1,0,0);
 MatrixSetIdentity(&m);
 MatrixRotate(&m,90,MatrixPrepend);
 AME(&m,0,1,-1,0,0,0);
 MatrixRotate(&m,270,MatrixPrepend);
 AME(&m,1,0,0,1,0,0);
 //
 MatrixSetValues(&m,10,20,30,40,50,60);
 MatrixRotateAt(&m,180,10,10,MatrixPrepend);
 AME(&m,-10,-20,-30,-40,850,1260);
 MatrixSetValues(&m,10,20,30,40,50,60);
 MatrixSetValues(&m2,10,20,30,40,50,60);
 MatrixMultiply(&m,&m,&m2);
 AME(&m,700,1000,1500,2200,2350,3460);
 MatrixSetValues(&m,10,20,30,40,50,60);
 MatrixSetValues(&m2,10,20,30,40,50,60);
 MatrixMultiply(&m,&m2,&m);
 AME(&m,700,1000,1500,2200,2350,3460);
 MatrixSetValues(&m,10,20,30,40,50,60);
 MatrixSetValues(&m2,10,20,30,40,50,60);
 ASSERT(MatrixEquals(&m,&m2));
 ASSERT(MatrixAreClose(&m,&m2));
 MatrixSetValue(&m,MatrixM11,10+1e-13);
 ASSERT(!MatrixEquals(&m,&m2));
 ASSERT(MatrixAreClose(&m,&m2));//may fail
 MatrixSetValue(&m,MatrixM11,20);
 ASSERT(!MatrixEquals(&m,&m2));
 ASSERT(!MatrixAreClose(&m,&m2));
 MatrixSetValues(&m,1,2,3,4,5,6);
 MatrixInvert(&m,&m);
 AME(&m,-2,1,1.5,-0.5,1,-2);
 MatrixSetValues(&m,1,0,0,1,8,8);
 MatrixInvert(&m,&m);
 AME(&m,1,0,0,1,-8,-8);
 MatrixSetIdentity(&m);
 MatrixInvert(&m,&m);
 ASSERT(MatrixIsIdentity(&m));
 ASSERT(MatrixIsInvertible(&m));
 MatrixSetValues(&m,10,20,30,40,50,60);
 MatrixScale(&m,2,4,MatrixPrepend);
 AME(&m,20,40,120,160,50,60);
 MatrixScale(&m,0.5,0.25,MatrixPrepend);
 AME(&m,10,20,30,40,50,60);
 MatrixSetValues(&m,5,10,15,20,25,30);
 MatrixScale(&m,2,4,MatrixAppend);
 AME(&m,10, 40, 30, 80, 50, 120);
 MatrixSetValues(&m,5,10,15,20,25,30);
 MatrixScaleAt(&m,2,4,100,100,MatrixAppend);
 AME(&m,10, 40, 30, 80, -50, -180);
 MatrixSetValues(&m,5,10,15,20,25,30);
 MatrixSkew(&m,45,180,MatrixAppend);
 AME(&m,15,10,35,20,55,30);
 MatrixSetValues(&m,5,10,15,20,25,30);
 MatrixRotateAt(&m,90,100,100,MatrixAppend);
 AME(&m,-10,5,-20,15,170,25);
 MatrixSetValues(&m,10,20,30,40,50,60);
 MatrixScale(&m,-2,-4,MatrixPrepend);
 AME(&m,-20,-40,-120,-160,50,60);
 MatrixSetValues(&m,10,20,30,40,50,60);
 MatrixShear(&m,2,4,MatrixPrepend);
 AME(&m,130,180,50,80,50,60);
 MatrixSetValues(&m,5,3,9,2,2,1);
 MatrixShear(&m,10,20,MatrixPrepend);
 AME(&m,185,43,59,32,2,1);
#define AP(x,a,b) \
  ASSERT(AreClose(points[x*2],a));\
  ASSERT(AreClose(points[x*2+1],b));
 {
 double points[]={2,4,4,8};
 MatrixSetValues(&m,2,4,6,8,10,12);
 MatrixTransformPoints(&m,points,2);
 AP(0,38,52);AP(1,66,92);
 }
 {
 double points[]={2,4,4,8};
 MatrixSetValues(&m,2,4,6,8,10,12);
 MatrixTransformVectors(&m,points,2);
 AP(0,28,40);AP(1,56,80);
 }
 {
 double rect[]={100,200,300,400};
 double points[]={10,20,30,40,50,60};
 MatrixParallelogram(&m,rect,points,MatrixSet);
 AME(&m,0.06666666,0.06666666,0.09999999,0.09999999,-16.6666679,-6.666667);
 }
 MatrixSetValues(&m,2,4,6,8,10,12);
 MatrixTranslate(&m,5,10,MatrixPrepend);
 AME(&m,2,4,6,8,80,112);
 DebugOut("End of tests");
 return 0;
}
*/
#endif

---