a+b as matrixes not seems ok, better c.sum(a, b)
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <string.h>
#include <float.h>
#include <iostream.h>
#define G(a,b) if(a)goto b
#define R return
#define P printf
#define M malloc
#define Fr free
#define F for
#define S sizeof
#define u8 unsigned char
#define u16 unsigned short
#define u32 unsigned
#define i32 int
#define d64 double
static u32 PBadAlloc;
};
template<class T> u32 permutazioni<T>::PBadAlloc=0;
ou32(u32 a){v=a;}
//operator u32() const{R v;}
ou32 operator=(u32 h){v=h; R *this;}
ou32 operator=( ou32 a){v=a.v; R *this;}
ou32 operator*=(ou32 a){v=omul(v, a.v); R *this;}
ou32 operator*=( u32 a){v=omul(v, a); R *this;}
ou32 operator+=(ou32 a){v=oadd(v, a.v); R *this;}
ou32 operator+=( u32 a){v=oadd(v, a ); R *this;}
ou32 operator-=(ou32 a){v=osub(v, a.v); R *this;}
ou32 operator-=( u32 a){v=osub(v, a ); R *this;}
ou32 operator/=(ou32 a){if(v==-1||a.v==-1) v=-1;
else v=v/a.v;
R *this;
}
ou32 operator/=(u32 a){if(v==-1||a==-1) v=-1;
else v=v/a;
R *this;
}
ou32 operator++(){v=(v==-1? -1: ++v); R *this;}
ou32 operator++(int){ou32 t=*this; if(v!=-1) ++v; R t;}
ou32 operator--(){v=(v==-1? -1: --v); R *this;}
ou32 operator--(int){ou32 t=*this; if(v!=-1) --v; R t;}
friend ou32 operator*(ou32 a, ou32 b){ou32 r; r.v=omul(a.v, b.v); R r;}
friend ou32 operator*(ou32 a, u32 b){ou32 r; r.v=omul(a.v, b ); R r;}
friend ou32 operator*(u32 a, ou32 b){ou32 r; r.v=omul(a , b.v); R r;}
i32 Mcon(u32** pm, u32 msize, ...);
u8* Aelm(u32* pm, ...);
template<class T> class Matrix{
public:
Matrix(u32 a, u32 b){Mcon(&mtr, S(T), a, b, -1);}
Matrix(u32 a, u32 b, u32 c){Mcon(&mtr, S(T), a, b, c, -1);}
Matrix(u32 a, u32 b, u32 c, u32 d){Mcon(&mtr, S(T), a, b, c, d, -1);}
Matrix(Matrix<T>& a)
{ if(a.mtr[0]==2)
R Mcon(&mtr, S(T), a.mtr[1], a.mtr[2], -1);
else if(a.mtr[0]==3)
R Mcon(&mtr, S(T), a.mtr[1], a.mtr[2], a.mtr[3], -1);
else if(a.mtr[0]==4)
R Mcon(&mtr, S(T), a.mtr[1], a.mtr[2], a.mtr[3], a.mtr[4], -1);
}
~Matrix(){Mcon(&mtr, 0);}
T& e(u32 i, u32 j) {R *(T*) Aelm(mtr, i,j);}
T& e(u32 i, u32 j, u32 k){R *(T*) Aelm(mtr, i,j,k); }
T& e(u32 i, u32 j, u32 k, u32 v){R *(T*) Aelm(mtr, i,j,k, v); }
i32 ov(u32 a){R 0;} // mai overflow
u32 ro(u32 a){R -1;}
i32 ov(i32 a){R 0;}
i32 ro(i32 a){R -1;}
i32 ov(d64 a){R _isnan(a)|| ! _finite(a) ? 1: 0; } // overflow con NAN
d64 ro(d64 a){d64 r=1.0E800; R r;} /* ritorna un elemento di d64 di errore*/
i32 ov(ou32 a){R a.v==-1? 1: 0; } // overflow con -1 range: 0..0FFFFFFFE
ou32 ro(ou32 a){ou32 r; r.v=-1; R r;} /* ritorna un elemento di ou32 di errore*/
i32 ov(oi32 a){R a.v==0x80000000? 1: 0; } // overflow con 0x80000000
/* ritorna un elemento di oi32 di errore*/
oi32 ro(oi32 a){oi32 r; r.v=0x80000000; R r;}
Matrix<T>& operator=(Matrix<T>& a)
{u32 i;
*Led(mtr)=0;
if(a.mtr== 0)
{
l0: *Led(mtr)=1; G(1,z); /* errore trovato e riportato */
}
G(*Led(a.mtr), l0);
G(*SZe(mtr)!=*SZe(a.mtr), l0);
if(a.mtr!=mtr)
{// deve essere delle stesse dimensioni e indici massimi
F(i=0;i<=mtr[0];++i)
G(a.mtr[i]!=mtr[i], l0);
memcpy(mtr, a.mtr, *SZ(a.mtr)+(2*a.mtr[0]+4)*sizeof(u32));
}
z:
R *this;
}
Matrix<T>& operator+=(Matrix<T>& a)
{u32 i;
T *p1,*p2;
u8 *pmax;
if(a.mtr== 0)
{
l0: *Led(mtr)=1; G(1,z); /* errore trovato e riportato */
}
G(*Led(mtr), z); G(*Led(a.mtr), l0);
// deve essere delle stesse dimensioni e indici massimi
F(i=0;i<=mtr[0];++i)
G(a.mtr[i]!=mtr[i], l0);
G(*SZe(mtr)!=*SZe(a.mtr), l0);
p1=(T*) Data(mtr); p2=(T*) Data(a.mtr);
i=*SZ(mtr); pmax=(u8*) p1; pmax+=i;
F(; p1<pmax ;++p1,++p2)
{*p1=*p1+*p2; G(ov(*p1),l0);}
z:
R *this;
}
Matrix<T>& operator-=(Matrix<T>& a)
{u32 i;
T *p1,*p2;
u8 *pmax;
if(a.mtr== 0)
{
l0: *Led(mtr)=1; G(1,z); /* errore trovato e riportato */
}
G(*Led(mtr), z); G(*Led(a.mtr), l0);
// deve essere delle stesse dimensioni e indici massimi
F(i=0;i<=mtr[0];++i)
G(a.mtr[i]!=mtr[i], l0);
*Led(mtr)=0;
G( *Led(a.mtr)||*Led(b.mtr), l0);
// la 1xn si puo' moltiplicare per nx1
// ovvero jxn si puo' moltiplicare per nxm
if(a.mtr[0]!=2||b.mtr[0]!=2||(*this).mtr[0]!=2)
{ l0: *Led(mtr)=1; l1: R -1; }
G( *SZe(mtr) != *SZe(a.mtr), l0);
G(a[1]!=b[0]||(*this)[0]!=a[0]||(*this)[1]!=a[1], l0);
if(this!= &a && this!= &b)
{i0=(*this)[0]; i1=(*this)[1];
F(i=0; i<i0; ++i)
{F(j=0; j<i1; ++j)
{F(z=0, s=0u; z<i1; ++z)
{a1=a.e(i,z); a2=b.e(z,j); s=s+(a1*a2);
G(ov(s),l0);
}
(*this).e(i,j)=s;
}
}
}
else {Matrix<T> t(mtr[1], mtr[2]);
t.mul(a,b); *this=t;
}
R 0;
}
// *this=a-b
i32 sub(Matrix<T>& a, Matrix<T>& b)
{u32 i;
T *p1,*p2, p3;
u8 *pmax;
*Led(mtr)=0;
if(a.mtr==0||b.mtr==0)
{
l0: *Led(mtr)=1; R -1;
}
G( *Led(a.mtr)||*Led(b.mtr), l0);
// deve essere delle stesse dimensioni e indici massimi
F(i=0;i<=mtr[0];++i)
G(a.mtr[i]!=mtr[i]||b.mtr[i]!=mtr[i], l0);
G(*SZe(mtr)!=*SZe(a.mtr), l0);
G(*SZe(mtr)!=*SZe(b.mtr), l0);
p1=(T*)Data(mtr); p2=(T*)Data(a.mtr); p3=(T*) Data(b.mtr);
i=*SZ(mtr); pmax=(u8*) p1; pmax+=i;
F(; p1<pmax; ++p1, ++p2)
{*p1=*p2-*p3; G(ov(*p1),l0);}
R 0;
}
// *this=a+b
i32 add(Matrix<T>& a, Matrix<T>& b)
{u32 i;
T *p1,*p2, p3;
u8 *pmax;
*Led(mtr)=0;
if(a.mtr==0||b.mtr==0)
{
l0: *Led(mtr)=1; R -1;
}
// deve essere delle stesse dimensioni e indici massimi
G( *Led(a.mtr)||*Led(b.mtr), l0);
F(i=0;i<=mtr[0];++i)
G(a.mtr[i]!=mtr[i]||b.mtr[i]!=mtr[i], l0);
G(*SZe(mtr)!=*SZe(a.mtr), l0);
G(*SZe(mtr)!=*SZe(b.mtr), l0);
p1=(T*)Data(mtr); p2=(T*)Data(a.mtr); p3=(T*) Data(b.mtr);
i=*SZ(mtr); pmax=(u8*) p1; pmax+=i;
F(; p1<pmax ;++p1,++p2)
{ *p1=*p2+*p3; G(ov(*p1),l0); }
R 0;
}
i32 transpose(Matrix<T>& a)
{u32 t, i, j, i0, i1;
*Led(mtr)=0;
if(a.mtr==0)
{
l0: *Led(mtr)=1; l1: R -1;
}
G(*Led(a.mtr), l0);
G(*SZe(mtr)!=*SZe(a.mtr), l0);
G(a.mtr[0]!=2||mtr[0]!=2, l0);
G(mtr[1]*mtr[2]!=a.mtr[1]*a.mtr[2], l0);
mtr[1]=a.mtr[2]; mtr[2]=a.mtr[1];
mtr[4]=a.mtr[4]; mtr[3]=a.mtr[1];
i0=mtr[1]; i1=mtr[2];
F(i=0; i<i0; ++i)
{F(j=0; j<i1; ++j)
(*this).e(i,j)=a.e(j,i);
}
R 0;
}
i32 transpose(void)
{u32 i, j, i0;
T t;
if(mtr[0]!=2)
{
l0: *Led(mtr)=1; l1: R -1;
}
G( *Led(mtr), l1 );
if(mtr[1]==mtr[2]) // NxN matrix
{i0=mtr[1];
F(i=0; i<i0; ++i)
{F(j=0; j<i; ++j)
{t=(*this).e(i,j);
(*this).e(i,j)=(*this).e(j,i);
(*this).e(j,i)=t;
}
}
}
else {Matrix<T> a(mtr[2], mtr[1]);
G(a.e(), l0);
G( a.transpose(*this)==-1 , l0);
i =mtr[1]; mtr[1]=mtr[2]; mtr[2]=i;
mtr[3]=i; // cambia gli indici per poter copiare
*this=a; // qui
}
R 0;
}
T det(void)
{u32 t, i1, k;
permutazioni<u32> pe(mtr[1]);
T rtn, p, s;
rtn=ro(rtn); // riempie rtn con valore errato
if(mtr[0]!=2u||pe.e())
{
l1: R rtn;
}
G( *Led(mtr), l1);
G(mtr[1]!=mtr[2], l1);
i1=mtr[1]; s=0;
la: F(k=0, p=1; k<i1; ++k)
p=p*(*this).e(k, pe[k]);
G(ov(p),l1);
if(pe.cnt%2!=0) {s=s-p; // cout<<"-"<<(i32)p;
}
else {s=s+p; // cout<<"+"<<(i32)p;
}
G(ov(s),l1);
G( (t=pe.nextP(1))==-1, l1);
G( t==1, la);
R s;
}
i32 e(void){R MtrBadAlloc;}
i32 ini(T* a)
{ T *p1, *p2;
u8 *pmax;
u32 i;
if(a==0) R -1;
*Led(mtr)=0;
p1=(T*) Data(mtr); p2=a;
i=*SZ(mtr); pmax=(u8*) p1; pmax+=i;
F(; p1<pmax ;++p1, ++p2)
*p1=*p2;
R 0;
}
/* ritorna il sup dell'indice consentito */
/* indici = 0, 1, 2 ,3 ,4 */
/* 2 i0 i1 => 2<=0 2<=1 2<=2* */
i32 operator[](i32 i)
{if(mtr==0||mtr[0]<=i||*Led(mtr)!=0) R 0;
R mtr[i+1];
}
friend Matrix<T> operator+(Matrix<T>& a, Matrix<T>& b)
{Matrix<T> r; r.add(a,b); R r;}
friend Matrix<T> operator-(Matrix<T>& a, Matrix<T>& b)
{Matrix<T> r; r.sub(a,b); R r;}
friend Matrix<T> operator*(Matrix<T>& a, Matrix<T>& b)
{Matrix<T> r; r.mul(a,b); R r;}
friend i32 operator==(Matrix<T>& a, Matrix<T>& b)
{T *p1, *p2;
u8 *pmax;
u32 i;
if( *Led(a.mtr)||*Led(b.mtr) )
R -1;
F(i=0; i<=a.mtr[0]; ++i)
if(a.mtr[i]!=b.mtr[i]) R 0;
p1=(T*)Data(a.mtr); p2=(T*) Data(b.mtr); i=*SZ(mtr);
if(i!=sizeof(T)) R -1;
pmax=(u8*) p1; pmax+=i;
F(; p1<pmax ;++p1,++p2)
if( *p1!=*p2 ) R 0;
R 1;
}
friend i32 operator!=(Matrix<T>& a, Matrix<T>& b){R !(a==b);}
friend ostream& operator<<(ostream& ostr, Matrix<T>& a)
{T *p1;
u8 *pmax;
u32 i, j;
ostr<<"\n";
if( *Led(a.mtr) ) {cout<<"[Me]\n"; R ostr;}
if( a.mtr[0]!=2 ) {cout<<"[M_Big]\n"; R ostr;}
p1=(T*)Data(a.mtr); i=*SZ(a.mtr);
pmax=(u8*) p1; pmax+=i; i=*SZe(a.mtr);
if(i>8||i==0) {cout<<"Out Error\n"; R ostr;}
F(j=0 ; p1<pmax ;++p1, ++j)
{ if(i==1) cout.width(3);
else if(i==2) cout.width(5);
else if(i==4) cout.width(10);
else if(i==8){cout.width(10); cout.precision(2);}
cout<< *p1;
if((j+1)%a.mtr[2]==0) cout<<"\n";
else cout<< " ";
}
R ostr;
}
friend istream& operator>>(istream& istr, Matrix<T>& a)
{T *p1;
u8 *pmax;
u32 i;
*Led(a.mtr)=0;
p1=(T*)Data(a.mtr); i=*SZ(a.mtr); pmax=(u8*) p1; pmax+=i;
F( ; p1<pmax ;++p1)
cin>>*p1;
R istr;
}
// for now only nxm matrix of u8,i8 u16,i16 u32,i32 d64[double]
void print(void)
{d64 dv;
u8 *p, *pmax;
u32 i,j,maxn, max, size, v;
u32 Index[256];
T a;
a=-1;
G(mtr, l1);
l0: R; // error
la: P("[Me]"); G(1,l0);
l1: G(mtr[0]!=2, l0);
size= *SZe(mtr);
G(size>8||size==0, l0);
G( *Led(mtr) , la);
i=0; p=Data(mtr); pmax=p+(size*mtr[1]*mtr[2]);
m0: G(size==1,l2);G(size==2,l3);G(size==4,l4);G(size==8,l11);G(1,l0);
l11: dv=*(d64*)p; P("%10.2f",dv); G(1,l5);
l2: v =*(u8*) p; G(a<0,l21); P("%3u ", (u32) v); G(1,l5);
l21: P("%3d ", (i32) v); G(1,l5);
l3: v =*(u16*)p; G(a<0,l31); P("%5u ", (u32) v); G(1,l5);
l31: P("%5d ", (i32) v); G(1,l5);
l4: v =*(u32*)p; G(a<0,l41); P("%10u ",(u32) v); G(1,l5);
l41: P("%10d ",(i32) v);
l5: G( (i+1)%mtr[2]!=0, l6); P("\n");
l6: p+=size; ++i; G(p<pmax, m0);
R;
}
void show(void)
{if(mtr)
{cout<<"n Index="<<mtr[0]<<"; i0="<<mtr[1]<<" ;i1="<<mtr[2];
cout<<"; w0="<<mtr[3]<<"; W1="<<mtr[4]<<"\n";
}
}
T A(u32 i, u32 j)
{T rtn;
u32 i1, i2, j1, j2;
rtn=ro(rtn); // riempie rtn con valore errato
if(*Led(mtr)) R rtn;
if(mtr[0]!=2|| mtr[1]!=mtr[2] ) R rtn;
if(i>=mtr[1]||j>=mtr[2]||mtr[1]<2) R rtn;
if(mtr[1]==2) // 2x2 1 elemento
{if(i==0 && j==0) R (*this).e(1,1);
if(i==0 && j==1) R (*this).e(1,0);
if(i==1 && j==0) R (*this).e(0,1);
if(i==1 && j==1) R (*this).e(0,0);
}
else {Matrix<T> a(mtr[1]-1, mtr[1]-1);
if(a.e()) R rtn;
j1=0; j2=0;
F(i1=0; i1<mtr[1]; ++i1)
F(i2=0; i2<mtr[1]; ++i2)
{if(i!=i1 && j!=i2)
{a.e(j1, j2++)=(*this).e(i1, i2);
if(j2==a.mtr[2]){j2=0; ++j1;}
}
}
rtn=a.det();
}
R rtn;
}
/* Return -1 if a has not inverse or error
Retrun 0 if make right calc with (*this)== a^(-1)
*/
i32 inverse(Matrix<T>& a)
{T adet, r;
u32 i,j;
*Led(mtr)=0;
if(mtr[0]!=2||a.mtr[0]!=2)
{
l0: *Led(mtr)=1;
R -1;
}
G( mtr[1]!=a.mtr[1]||mtr[2]!=a.mtr[2]||mtr[1]!=mtr[2] , l0);
adet=2; adet/=2;
G(adet==0, l0); /* il tipo non ammette inversa */
adet=a.det();
G(adet==0u, l0);
F(i=0; i<mtr[1]; ++i)
F(j=0; j<mtr[1]; ++j)
{r=((i+j)%2==0? +1: -1 )* a.A(j,i)/adet;
G(ov(r), l0);
(*this).e(i,j)=r;
}
R 0;
}
u32 *mtr;
static u32 MtrBadAlloc;
};
template<class T> u32 Matrix<T>::MtrBadAlloc=0;
/*
reserve mem space to matrix of element msize for dimension
i0xi1x...xin taken as input
where i0, i1 etc are max+1 index and n<=255
the last max+1 index is -1 that end the input
Return -1 for error
Return 0 all ok
NB
**The first function to call on pm [construct matrix] **
**the last function to call on pm [destruct matrix] **
*/
i32 Mcon(u32** pm, u32 msize, ...)
{va_list ap;
u32 v, i, s, j;
u32 arr[256];
G(pm, l1);
l0: Matrix<u32>::MtrBadAlloc=1; R -1; // for error
u32 f(void)
{Matrix<ou32> m1(2,2),m2(2,3),m6(2,2),m7(2,2),m8(1,2),m9(1,2);
Matrix<d64> m3(2,2), md1(2,2), md2(2,2), mdd(2,2), md3(3,3), md4(3,3);
Matrix<ou32> m10(2,2);
Matrix<u8> m4(1,7), m5(7,1);
Matrix<oi32> m12(4,4), m13(3,3);
Matrix<oi32> m20(4,2), m21(2,4), m23(2,4), m24(2,4);
u32 d, i, j;
ou32 dd;
i32 vm[] ={ 2,9,9,4, 2,-3,12,8, 4,8,3,-5, 1,2,6,4 }, r;
i32 vm1[]={ 1,2,3, 4,5,6, 7,8,9};
/* if some Matrix<T> object is not allocated return 0 */
if(m1.e()) R 0;
F(i=0; i<m2[0]; ++i)
F(j=0; j<m2[1]; ++j)
m2.e(i,j)=10*i+j;
F(i=0; i<m4[0]; ++i)
F(j=0; j<m4[1]; ++j)
m4.e(i,j)=10*i+j;
F(i=0; i<m5[0]; ++i)
F(j=0; j<m5[1]; ++j)
m5.e(i,j)=10*i+j;
m1.e(0,0)=1; m1.e(0,1)=2;
m1.e(1,0)=3; m1.e(1,1)=4;
dd = m1.det();
m6.inverse(m1);
m12.ini( (oi32*) vm ); /* can be wrong etc */
m13.ini( (oi32*) vm1);
m6.e(0,0)=5; m6.e(0,1)=6;
m6.e(1,0)=7; m6.e(1,1)=8;
m7.mul(m1, m6);
F(i=0; i<m3[0]; ++i)
F(j=0; j<m3[1]; ++j)
m3.e(i,j)=10.23*i+j;
m8.e(0,0)=1; m8.e(0,1)=2;
m9.mul(m8, m1);
P("M1:\n"); m1.print(); P("\n\n");
P("m1.det()=%d\n", d);
P("M2:\n"); m2.print(); P("\n\n");
P("M3:\n"); m3.print(); P("\n\n");
P("M4:\n"); m4.print(); P("\n\n");
m4.show();
P("M5:\n"); m5.print(); P("\n\n");
m5.show();
P("M6:\n"); m6.print(); P("\n\n");
P("M7=m1*m6:\n"); m7.print(); P("\n\n");
P("M8\n"); m8.print(); P("\n\n");
P("M9=m8*m1\n"); m9.print(); P("\n\n");
P("M12\n");m12.print(); P("\n\n");
r=m12.det();
P("DET(M12)=%d\n", r);
F(i=0; i<m12.mtr[1]; ++i)
F(j=0; j<m12.mtr[2]; ++j)
cout<< "A("<<i<<", "<< j<< ")="<< m12.A(i,j) << " ";
m7=m1;
P("M7: m7=m1\n"); m7.print(); P("\n\n");
m1*=(ou32)10u;
P("M1=M1*10:\n"); m1.print(); P("\n\n");
m1+=m1;
P("M1+=M1:\n"); m1.print(); P("\n\n");
P("M13:\n"); m13.print(); P("\n\n");
r=m13.det();
P("DET(M13)=%d\n", r);
m1.transpose();
P("M1 traspose:\n"); m1.print(); P("\n\n");
m20.ini( (oi32*) vm1); /* can be wrong etc */
m21.ini( (oi32*) vm1);
P("M20:\n"); m20.print(); P("\n"); m20.show(); P("\n\n");
P("M21:\n"); m21.print(); P("\n"); m21.show(); P("\n\n");
m23.transpose(m20);
P("M23:\n"); m23.print(); P("\n"); m23.show(); P("\n\n");
m24.transpose(m21);
P("M24:\n"); m24.print(); P("\n"); m24.show(); P("\n\n");
m20.transpose();
P("M20:\n"); m20.print(); P("\n"); m20.show(); P("\n\n");
m21.transpose();
P("M21:\n"); m21.print(); P("\n"); m21.show(); P("\n\n");
cout<<"M1 M2="<< m1 << m21;
//cin>>m1; cout<<m1;
m6.inverse(m1);
cout<<"M1-1==\n"<<m6<<"\n\n";
F(i=0; i<md1.mtr[1]; ++i)
F(j=0; j<md1.mtr[2]; ++j)
md1.e(i, j)=10.0*i+j;
cout<<"Md1:"<< md1;
md2.inverse(md1);
cout<<"Md1-1"<<md2;
md1.mul(md1, md2);
cout<<"I="<< md1;
d=70;
F(i=0; i<md3.mtr[1]; ++i)
F(j=0; j<md3.mtr[2]; ++j)
md3.e(i, j)=d++;
++md3.e(1,1);
cout<<"Md3:"<< md3;
cout<<"|Md3|="<< md3.det()<<" ";
md4.inverse(md3);
cout<<"Md3-1"<<md4;
md3.mul(md3, md4);
cout<<"I="<< md3;