Column slices

[luser droog]

I've been thinking further about the arbitrary-dimensional array
structure (and am still stunned by the lackluster response so far),
and I think a suitable coup-de-grace for demonstrating the value
of indirection is to tease-out yet another data member for this
structure. This ties everything back to the 1962 APL book, the
chapter on Representation, which has kept me up so many, many nights.

A dynamic array object is defined by its rank, list of dimension sizes,
a pointer to the start of the data, AND a weighting vector which for
simple uses merely caches the multiplications from the array-access
formula but it affords us much greater flexibility as (I hope) we
shall see...

typedef struct arr {
int rank;
int *dims;
int *weight;
int *data;
} *parr;

All of the pointer members can be assigned locations within the
same allocated block of memory as the struct itself (which we will
call the header). But by replacing the earlier use of offsets
and struct-hackery, the implementation of algorithms can be made
independent of that actual memory layout within (or without) the
block.

The basic memory layout for a self-contained array object is

rank dims weight data
dims[0] dims[1] ... dims[rank-1]
weight[0] weight[1] ... weight[rank-1]
data[0] data[1] ... data[ product(dims)-1 ]

An indirect array sharing data (whole array or 1 or more row-slices)
will have the following memory layout

rank dims weight data
dims[0] dims[1] ... dims[rank-1]
weight[0] weight[1] ... weight[rank-1]
//no data! it's somewhere else!

And an indirect array containing an orthogonal slice along
another axis will have the same layout as a basic indirect array,
but with dims and weight suitably modified.

The access formula for an element with indices (i0 i1 ... iN)
is

pa->data[ i0*pa->weight[0] + i1*pa->weight[1] + ...
+ iN*pa->weight[N] ]

, assuming each index i[j] is between [ 0 ... dims[j] ).

Nothing implemented yet.

[asetof...@gmail.com]

Is weight[i] == dim[i]*sizeof(matrix element)?

[Nobody]

On Tue, 05 May 2015 23:18:12 -0700, luser droog wrote:

> A dynamic array object is defined by its rank, list of dimension sizes, a
> pointer to the start of the data, AND a weighting vector which for simple
> uses merely caches the multiplications from the array-access formula but
> it affords us much greater flexibility as (I hope) we shall see...

This is essentially NumPy's ndarray (N-dimensional array) class:

http://docs.scipy.org/doc/numpy/reference/arrays.ndarray.html

That may be worth looking at if you're going to reinvent that particular
wheel (the concept certainly pre-dates NumPy, but it's well documented,
maintained, and doesn't require a fancy keyboard).

[Bartc]

On 06/05/2015 07:18, luser droog wrote:

> A dynamic array object is defined by its rank, list of dimension sizes,
> a pointer to the start of the data, AND a weighting vector which for
> simple uses merely caches the multiplications from the array-access
> formula but it affords us much greater flexibility as (I hope) we
> shall see...
>
> typedef struct arr {
> int rank;
> int *dims;
> int *weight;
> int *data;
> } *parr;
>
> All of the pointer members can be assigned locations within the
> same allocated block of memory as the struct itself (which we will
> call the header). But by replacing the earlier use of offsets
> and struct-hackery, the implementation of algorithms can be made
> independent of that actual memory layout within (or without) the
> block.
>
> The basic memory layout for a self-contained array object is
>
> rank dims weight data
> dims[0] dims[1] ... dims[rank-1]
> weight[0] weight[1] ... weight[rank-1]
> data[0] data[1] ... data[ product(dims)-1 ]

So for a 3x4x5 array it would be:

3 4 5
W0 W1 W2
D0 D1 ..... D59 ?

If the last index (the 0 to 4 one) varies more quickly, the data
consists of 3x4 rows of 5, or 3 blocks of (4 rows of 5). And the data
will be more like:

D(0,0,0) D(0,0,1),... D(0,0,4) D(0,1,0) .... D(2,3,3) D(2,3,4)

Presumably the W vector might be W0=12 W1=4 W2=1, to help locate
D[i][j][k], with the last W not really needed.

But where does the column slice come into it? To access any slice other
than a row, would need the ability to deal with an 'array' that has a
stride, or offset between elements, the stride being W0, W1 (or W2 for a
row).

>
> An indirect array sharing data (whole array or 1 or more row-slices)
> will have the following memory layout
>
> rank dims weight data
> dims[0] dims[1] ... dims[rank-1]
> weight[0] weight[1] ... weight[rank-1]
> //no data! it's somewhere else!
>
> And an indirect array containing an orthogonal slice along
> another axis will have the same layout as a basic indirect array,
> but with dims and weight suitably modified.

Don't forget 2D, 3D and N-dimensional slices, where the slice (for a 3D
array) might be rectangular 2D array (in, I think, 3 possible
orientations), or a cuboid, which is a 3D subarray (and might also come
in different orientations, although you can't tell apart from the
indexing order).

> Nothing implemented yet.

> typedef struct arr {
> int rank;
> int *dims;
> int *weight;
> int *data;
> } *parr;

Does it need a different struct depending on the array element type?
This also seems inefficient for small dimensions, for 2D for example
which perhaps is as high as is usually needed.

And for 1D, although I guess a different method is used for those,
unless it is possible for the number of dimensions to be variable? Ie.
be itself a (1D) vector, requiring a similarly sized vector of indices.
(I suppose an N-dimensional array of dimensions is a possibility but I
can't get my head around that!)

(I guess these arrays form part of a different language (APL?), and C is
used to implement it. For such purposes, I don't bother with
N-dimensional arrays at all, just have arrays of arrays to any degree.
It's probably not quite as efficient, but it's simpler. And there are no
multiplications involved.)

--
Bartc

[luser droog]

On Wednesday, May 6, 2015 at 2:26:12 AM UTC-5, asetof...@gmail.com wrote:
> Is weight[i] == dim[i]*sizeof(matrix element)?

It's the product of all "lower" dimensions.
weight[i] = product(dims[i+1 .. rank-1]

So for dims = (3, 4, 5)
weight = (20, 5, 1) = (4*5, 5, 1)

or for dims = (7, 6, 5, 4, 3, 2)
weight = (720, 120, 24, 6, 2, 1)

[Bartc]

On 06/05/2015 10:30, Bartc wrote:
> So for a 3x4x5 array it would be:

> Presumably the W vector might be W0=12 W1=4 W2=1, to help locate

Or W0=20 as you demonstrate above. (The boxes in the diagram I created
were too badly drawn to see properly what size they were!)

--
Bartc

[luser droog]

I get confused when I try to do it directly, but you can get column
slices by taking a transpose and then taking row slices.

So if you have (d,w) = ((3,4,5),(20,5,1), then (d,w)' = ((5,3,4),(1,20,5))
and (d,w)'[0;*;*] = ((1,3,4)(1,20,5)). The transpose of this resulting
array is the [*;*;0] column slice of the original (d,w).

The weighting vector AFAIK is the sine qua non of taking column slices.
And yes, the last weight will usually be 1, unless you're working with
a transpose.

> >
> > An indirect array sharing data (whole array or 1 or more row-slices)
> > will have the following memory layout
> >
> > rank dims weight data
> > dims[0] dims[1] ... dims[rank-1]
> > weight[0] weight[1] ... weight[rank-1]
> > //no data! it's somewhere else!
> >
> > And an indirect array containing an orthogonal slice along
> > another axis will have the same layout as a basic indirect array,
> > but with dims and weight suitably modified.
>
> Don't forget 2D, 3D and N-dimensional slices, where the slice (for a 3D
> array) might be rectangular 2D array (in, I think, 3 possible
> orientations), or a cuboid, which is a 3D subarray (and might also come
> in different orientations, although you can't tell apart from the
> indexing order).

I think all those are possible, so long as the data can be accessed "linearly"
in the algebraic sense.

> > Nothing implemented yet.
>
> > typedef struct arr {
> > int rank;
> > int *dims;
> > int *weight;
> > int *data;
> > } *parr;
>
> Does it need a different struct depending on the array element type?
> This also seems inefficient for small dimensions, for 2D for example
> which perhaps is as high as is usually needed.

That gave me a lot of grief with inca2 which has 3 basic datatypes:
char, int, and double. It was indeed a royal pain to deal with the
sizes all over the place. For the new rewrite, I had a different
idea to encode everything into int-size handles. This is described
in this post in comp.lang.apl:
https://groups.google.com/d/topic/comp.lang.apl/KSwo40-eLKI/discussion

> And for 1D, although I guess a different method is used for those,
> unless it is possible for the number of dimensions to be variable? Ie.
> be itself a (1D) vector, requiring a similarly sized vector of indices.
> (I suppose an N-dimensional array of dimensions is a possibility but I
> can't get my head around that!)

Yes, the number of dimensions is variable and is stored in the .rank
struct member. The dimension list is 1D, a vector.

> (I guess these arrays form part of a different language (APL?), and C is
> used to implement it. For such purposes, I don't bother with
> N-dimensional arrays at all, just have arrays of arrays to any degree.
> It's probably not quite as efficient, but it's simpler. And there are no
> multiplications involved.)
>

Yes, this all comes directly from my experience trying to implement a
(better) APL interpreter in C. But I hope the ideas have more general
applicability. And perhaps I can build a library or at least a demo
of some linear algebra something or other that may help somebody
somewhere someday.

[asetof...@gmail.com]

Thank you, clear

[asetof...@gmail.com]

On Thursday, May 7, 2015 at 8:02:35 PM UTC+2, asetof...@gmail.com wrote:
> Thank you, clear


Here my start in implementing matrices ...
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <string.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

/*
I[0],I[1],I[2],I[3],I[4],sizeElement, sz, led, [matrix data]
0 1 2 3 4 5 6 7 8
1 2 3 4 5 6 7 8 9
I[I[0]+1]=I[3]
I[I[0]+I[0]]=I[4]

I[0],I[1],..,I[p[0]], sizeElement, sz, led, [matrix data]
0 1 I[0] 2*I[0]+1 2*I[0]+2 2*I[0]+3 2*I[0]+4 indirizzo
1 2 I[0]+1 2*I[0]+2 2*I[0]+3 2*I[0]+4 2*I[0]+5 nelement
I[0]: is the number of index as 32 bit
I[1]: is the first max index as 32 bit
....
I[I[0]]: is the last max index
I[I[0]+1]: is w[0]
...
I[2*I[0]]: is w[I[0]]
I[2*I[0]+1]: is the size of element as 32 bit 5
I[2*I[0]+2]: is sz size of "data" area 6
I[2*I[0]+3]: is led: 0 matrix is ok, 1 overflow or errors in matrix 7
sono 2*I[0]+4 elementi [8 elementi]
I[2*I[0]+4]: is matrix "data" area 8
*/

// access to led as u32*
#define Led(a) ((a)+*(a)*2+3)
// access to size of matrix data as u32*
#define SZ(a) ((a)+*(a)*2+2)
// access to size element of matrix as u32*
#define SZe(a) ((a)+*(a)*2+1)
// access to data as u8*
#define Data(a) (u8*) ((a)+*(a)*2+4)

/*
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: R -1; // for error
n0: va_end(ap); *pm=0; G(1,l0);
l1: G(msize, l3);Fr(*pm); *pm=0;/* free the Matrix */
l2: R 0; // for all ok
y0: va_end(ap); G(1,l2);
l3: va_start(ap, msize); G(msize>0xFFF, n0); G(msize==0, n0); i=0;
l4: v=va_arg(ap, u32); //v the first index;if i==0 means no element
G(v>0xFFFF && v!= -1, n0); G(v!=-1, la); G(i==0, n0);
s=1; j=1;
l5: s*=arr[j-1]; G(s>0xFFFFF, n0); ++j; G( j<=i, l5);
s*=msize; G(s>0xFFFFF, n0);
*pm=(u32*)M((2*i+4)*sizeof(u32)+s); G(*pm==0, n0);
(*pm)[0]=i; *SZe(*pm)=msize; *SZ(*pm)=s; *Led(*pm)=0;
j=1; s=1;
l6: (*pm)[j]=arr[j-1]; (*pm)[2*i+1-j]=s;
s*=arr[i-j]; ++j; G(j<=i, l6);
G(1,y0);
la: G(i>255, n0); arr[i++]=v; G(1,l4);
}

// return one pointer to element
// or 0 if find some error
/*
2x2 y*n +x (y, x)e{0..m -1}x{0..n -1}
2 elm(i0,i1)=i0*M1+i1 (i0,i1)e{0..M0-1}x{0..M1-1} Max=(M0-1)*M1+M1-1=M0*M1-1
se M0=7 M1=1 7x1 (1,1)
elm(i0,i1)=i0+i1 (i0,i1)e{0..6}x{0..0} Max=(7-1)*1+1-1=6
se M0=1 M1=7 1x7 (7,1)
elm(i0,i1)=i0*7+i1 (i0,i1)e{0..0}x{0..6} Max=(0)*7+7-1=6

3 elm(i0,i1,i2)=i0*M1*M2+i1*M1+i2 (i0,i1,i2)e{0..M0-1}x{0..M1-1}x{0..M2-1}
Max=(M0-1)*M1*M2+(M1-1)*M2+(M2-1)=M0*M1*M2-M1*M2+M1*M2-M2+M2-1=M0*M1*M2-1
*/
u8* Aelm(u32* pm, ...)
{va_list ap;
u8 *r;
u32 maxn, size, v, s, i;

G(pm,l1);
l0: R 0;
n0: va_end(ap); G(1,l0);
l1: maxn=pm[0]; size= *SZe(pm); G(size>0xFFFF, l0);
// 0 1
G(maxn>=256||maxn==0,l0); //for 2x2 y*n+x (y,x)e[0..m-1]x[0..n-1]
va_start(ap, pm); // gets all index i0 i1 maxn=2
i=1; s=0; /* I[0]+1..2*I[0] */
l3: v=va_arg(ap,u32); G(v>=pm[i],n0); s+=v*pm[maxn+i]; ++i; G(i<=maxn, l3);
v=s*size; r=Data(pm); r+=v;
va_end(ap);
R r;
}

// for now only nxm matrix of u8 u16 u32 d64[double]
void PrntMtrx(u32* m)
{va_list ap;
d64 dv;
u8 *p, *pmax;
u32 i,j,maxn, max, size, v;
u32 Index[256];

G(m, l1);
l0: R; // error
la: P("[Me]"); G(1,l0);
l1: G(m[0]>2, l0);
size= *SZe(m);
G(size>8||size==0, l0);
G( *Led(m) , la);
i=0; p=Data(m); pmax=p+(size*m[1]*m[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; P("%3u ", v); G(1,l5);
l3: v =*(u16*)p; P("%5u ", v); G(1,l5);
l4: v =*(u32*)p; P("%10u ", v);
l5: G( (i+1)%m[2]!=0, l6); P("\n");
l6: p+=size; ++i; G(p<pmax, m0);
R;
}

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(){Mcon(&mtr, 0);}

T& elm(u32 i, u32 j){R *(T*) Aelm(mtr, i,j);}
T& elm(u32 i, u32 j, u32 k){R *(T*) Aelm(mtr, i,j,k); }
T& elm(u32 i, u32 j, u32 k, u32 v){R *(T*) Aelm(mtr, i,j,k, v); }

Matrix<T>& operator=(Matrix<T>& a)
{u32 i;
if(a.mtr== 0)
{
l0: *Led(mtr)=1; G(1,z); /* errore trovato e riportato */
}
if(a.mtr!=mtr)
{// deve essere delle stesse dimensioni e indici massimi
G(a.mtr[0]!=mtr[0], l0);
F(i=1;i<=mtr[0];++i)
G(a.mtr[i]!=mtr[i], l0);
G(*Led(a.mtr)!=0, l0); G(*SZe(mtr)!=*SZe(a.mtr), 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 */
}
// deve essere delle stesse dimensioni e indici massimi
G(a.mtr[0]!=mtr[0], l0);
F(i=1;i<=mtr[0];++i)
G(a.mtr[i]!=mtr[i], l0);
G(*Led(a.mtr)!=0, 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;
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 */
}
// deve essere delle stesse dimensioni e indici massimi
G(a.mtr[0]!=mtr[0], l0);
F(i=1;i<=mtr[0];++i)
G(a.mtr[i]!=mtr[i], l0);
G(*Led(a.mtr)!=0, 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;
z:
R *this;
}

Matrix<T>& operator*=(T a)
{u32 i;
T *p1;
u8 *pmax;

p1=(T*) Data(mtr);
i=*SZ(mtr); pmax=(u8*) p1; pmax+=i;
F(; p1<pmax ;++p1) *p1 = *p1 * a;
R *this;
}

/* ritorna il sup dell'indice consentito */
/* indici = 0, 1, 2 ,3 ,4 */
u32 operator[](u32 i){if(mtr==0||mtr[0]<=i) R 0; R mtr[i+1];}
int print(void){PrntMtrx(mtr);}

u32 *mtr;
};

int main(void)
{Matrix<u32> m1(2,2), m2(2,3), m7(2,2);
Matrix<d64> m3(2,2);
Matrix<u8> m4(1,7), m5(7,1);
u32 i,j;

F(i=0; i<m2[0]; ++i)
F(j=0; j<m2[1]; ++j)
m2.elm(i,j)=10*i+j;

F(i=0; i<m4[0]; ++i)
F(j=0; j<m4[1]; ++j)
m4.elm(i,j)=10*i+j;

F(i=0; i<m5[0]; ++i)
F(j=0; j<m5[1]; ++j)
m5.elm(i,j)=10*i+j;

m1.elm(0,0)=123; m1.elm(0,1)=124;
m1.elm(1,0)=125; m1.elm(1,1)=126;

F(i=0; i<m3[0]; ++i)
F(j=0; j<m3[1]; ++j)
m3.elm(i,j)=10.23*i+j;

P("M1:\n"); m1.print(); P("\n\n");
P("M2:\n"); m2.print(); P("\n\n");
P("M3:\n"); m3.print(); P("\n\n");
P("M4:\n"); m4.print(); P("\n\n");
P("M5:\n"); m5.print(); P("\n\n");

m7=m1;
P("M7:\n"); m7.print(); P("\n\n");
m1*=10;
P("M1=M1*10:\n"); m1.print(); P("\n\n");
m1+=m1;
P("M1+=M1:\n"); m1.print(); P("\n\n");

R 0;
}

----------------------
M1:
123 124
125 126


M2:
0 1 2
10 11 12


M3:
0.00 1.00
10.23 11.23


M4:
0 1 2 3 4 5 6


M5:
0
10
20
30
40
50
60


M7:
123 124
125 126


M1=M1*10:
1230 1240
1250 1260


M1+=M1:
2460 2480
2500 2520

[asetof...@gmail.com]

"The dimension list is 1D, a vector"

What about
1D not exist in matrix?

I in my code not allow that
it does not exist.

Exist nx1
As column matrix for example for 3x1
1
2
3
Or 1xn matrix example 1x3
1 2 3
For both program use 2 index not only one
so is what you call 'rank' 2?

I forget all math on matrix
and follow the computer...

[luser droog]

A basic implementation. Insufficient error checking.

//arrind.c - Arrays, indirect.

#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "ppnarg.h"

typedef struct arr {
int rank;
int *dims;
int *weight;
int *data;

} *arr;

/* multiply together rank integers in dims array
*/
int productdims(int rank, int *dims){
int i,z=1;
for(i=0; i<rank; i++)
z *= dims[i];
return z;
}

/* create a new array with specified dimensions */
#define array(...) (makearr)(PP_NARG(__VA_ARGS__),__VA_ARGS__)

/* create a new array with specified rank and dimensions
*/
arr (makearr)(int rank, ...){
va_list ap;
//int *dims=calloc(rank,sizeof(int));
int dims[rank];
int datasz;
int i;
int x;
arr z;

va_start(ap,rank);
for (i=0; i<rank; i++){
dims[i]=va_arg(ap,int);
}
va_end(ap);

datasz=productdims(rank,dims);
z=malloc(sizeof(struct arr)
+ (rank+rank+datasz)*sizeof(int));

z->rank = rank;
z->dims = (int*)(((char*)z) + sizeof(struct arr));
z->weight = z->dims + rank;
z->data = z->weight + rank;
memmove(z->dims,dims,rank*sizeof(int));
for(x=1, i=rank-1; i>=0; i--){
z->weight[i] = x;
x *= z->dims[i];
}

//free(dims);
return z;
}

/* create a new array which shares the data of an existing array
*/
arr clone(arr a){
arr z=malloc(sizeof(struct arr)
+ (a->rank+a->rank)*sizeof(int));
z->rank = a->rank;
z->dims = (int*)(((char*)z) + sizeof(struct arr));
z->weight = z->dims + z->rank;
z->data = a->data;
memmove(z->dims,a->dims,z->rank*sizeof(int));
memmove(z->weight,a->weight,z->rank*sizeof(int));
return z;
}

/* exchange the leftmost two dimensions (only two in 2D)
*/
void transpose2(arr a){
int t;
//if (a->rank != 2) error();
t = a->dims[0]; a->dims[0] = a->dims[1]; a->dims[1] = t;
t = a->weight[0]; a->weight[0] = a->weight[1]; a->weight[1] = t;
}

/* take a (row) slice (in 2D) */
arr slice(arr a, int i){
int rank = a->rank-1;
arr z=malloc(sizeof(struct arr)
+ (rank+rank)*sizeof(int));
z->rank = rank;
z->dims = (int*)(((char*)z) + sizeof(struct arr));
z->weight = z->dims + z->rank;
memmove(z->dims,a->dims+1,z->rank*sizeof(int));
memmove(z->weight,a->weight+1,z->rank*sizeof(int));
z->data = a->data + i*a->weight[0];
return z;
}

/* access element of a indexed by int[] */
int *elema(arr a, int *ind){
int idx = 0;
int i;
for (i=0; i<a->rank; i++){
idx += ind[i] * a->weight[i];
}
return a->data + idx;
}

/* access element of a indexed by va_list */
int *elemv(arr a, va_list ap){
int idx = 0;
int i;
for(i=0; i<a->rank; i++){
int ind;
ind = va_arg(ap,int);
//if (ind > a->dims[i]) error();
idx += ind * a->weight[i];
}
return a->data + idx;
}

/* access element of a indexed by integer arguments */
int *elem(arr a, ...){
va_list ap;
int *z;

va_start(ap,a);
z = elemv(a,ap);
va_end(ap);

return z;
}

/* create a (contiguous) copy of a (not necessarily contiguous)
existing array
*/
arr copy(arr a){
int datasz = productdims(a->rank,a->dims);
arr z=malloc(sizeof(struct arr)
+ (a->rank+a->rank+datasz)*sizeof(int));
int i;
int x;
int ind[a->rank];

z->rank = a->rank;
z->dims = (int*)(((char*)z) + sizeof(struct arr));
z->weight = z->dims + z->rank;
z->data = z->weight + z->rank;
memmove(z->dims,a->dims,z->rank*sizeof(int));
for(x=1, i=z->rank-1; i>=0; i--){
z->weight[i] = x;
x *= z->dims[i];
}

for (i=0;i<datasz;i++){
int j;
int idx = i;
for (j=z->rank-1; j>=0; j--) { //generate index
ind[j] = idx % z->dims[j];
idx /= z->dims[j];
}
z->data[i] = *elema(a,ind);
}

return z;
}

#define OPERATORS(_) \
/* f F id */ \
_('+',+,0) \
_('*',*,1) \
_('=',==,1) \
/**/

/* perform binary operation F upon corresponding elements of vectors X and Y */
#define binop(X,F,Y) (binop)(X,*#F,Y)
#define BINOP(f,F,id) case f: *elem(z,i) = *elem(x,i) F *elem(y,i); break;
arr (binop)(arr x, char f, arr y){
arr z=copy(x);
int n=x->dims[0];
int i;
for (i=0; i<n; i++){
switch(f){
OPERATORS(BINOP)
}
}
return z;
}
#undef BINOP

/* perform binary operation F upon adjacent elements of vector X, right to left,
reducing vector to a single value */
#define reduce(F,X) (reduce)(*#F,X)
#define REDID(f,F,id) case f: x = id; break;
#define REDOP(f,F,id) case f: x = *elem(a,i) F x; break;
int (reduce)(char f, arr a){
int n = a->dims[0];
int x;
int i;
switch(f){
OPERATORS(REDID)
}
if (n) {
x=*elem(a,n-1);
for (i=n-2;i>=0;i--){
switch(f){
OPERATORS(REDOP)
}
}
}
return x;
}
#undef REDID
#undef REDOP

/* perform a (2D) matrix multiplication upon x and y
using operations F and G.
Z[i,j] = F/ X[i,*] G Y'[j,*]
*/
#define matmul(X,F,G,Y) (matmul)(X,*#F,*#G,Y)
arr (matmul)(arr x, char f, char g, arr y){
int i,j;
arr z=array(x->dims[0],y->dims[1]);

y=clone(y);
transpose2(y);
for (i=0; i<x->dims[0]; i++){
for (j=0; j<y->dims[0]; j++){
arr xs = slice(x,i);
arr ys = slice(y,j);
*elem(z,i,j)=(reduce)(f,(binop)(xs,g,ys));
free(xs);
free(ys);
}
}

free(y);
return z;
}

int main(){

#if 0
/* testing basic functionality and copies, transposes, and slices */
int loop;
for (loop = 0;
loop < 1/*00000*/;
loop++)
{

{
int i,n=12;
arr a=array(n);

for (i=0;i<n;i++)
*elem(a,i) = n-i;
for (i=0;i<n;i++,printf(" "))
printf("%2d",*elem(a,i));
printf("\n\n");

free(a);
}

{
int i,j,n=6;
arr a=array(n,n);
arr b;

for (i=0;i<n;i++)
for (j=0;j<n;j++)
*elem(a,i,j) = n*n - (i*n+j);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf(" "))
printf("%2d",*elem(a,i,j));
printf("\n");

b=clone(a);
transpose2(b);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf(" "))
printf("%2d",*elem(b,i,j));
printf("\n");

free(b);
free(a);
}

{
int i,j,k,n=3;
arr a=array(n,n,n);
arr b;
arr c;
arr d;
for (i=0;i<n;i++)
for (j=0;j<n;j++)
for (k=0;k<n;k++)
*elem(a,i,j,k) = n*n*n - ((i*n+j)*n+k);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf("\n"))
for (k=0;k<n;k++,printf(" "))
printf("%2d",*elem(a,i,j,k));
printf("\n");

/*
b=slice(a,0);
transpose2(b);
for (i=0;i<n;i++)
for (j=0;j<n;j++)
*elem(b,i,j) = n*n - (i*n+j);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf("\n"))
for (k=0;k<n;k++,printf(" "))
printf("%2d",*elem(a,i,j,k));
printf("\n");
free(b);
*/

b=clone(a);
transpose2(b);
c=slice(b,1);
transpose2(c);
for (i=0;i<n;i++)
for (j=0;j<n;j++)
*elem(c,i,j) = n*n - (i*n+j);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf("\n"))
for (k=0;k<n;k++,printf(" "))
printf("%2d",*elem(b,i,j,k));
printf("\n");

//d = copy(a);
d = copy(b);
for (i=0; i<n*n*n; i++){
printf("%2d", d->data[i]);
printf(" ");
if (((i+1)%n)==0) printf("\n");
if (((i+1)%(n*n))==0) printf("\n");
}

free(a);
free(b);
free(c);
free(d);
}
}
#endif

{ /* testing reduce() */
int i,n=3;
arr a=array(n);
arr b;
for (i=0; i<n; i++)
*elem(a,i) = i+1;
printf("%2d\n", reduce(*,a));

free(a);

n=6;
a=array(n);
for (i=0; i<n; i++)
*elem(a,i) = 5;
b=binop(a,=,a);
for (i=0; i<n; i++,printf(" "))
printf("%2d", *elem(b,i));
printf("\n%2d\n", reduce(=,b));
free(a);
free(b);
}

{ /* testing matmul() */
int i,j,n=3;
arr a=array(n,n);
arr b;

for (i=0;i<n;i++)
for (j=0;j<n;j++)
*elem(a,i,j) = ((i*n)+j)+1;
b = matmul(a,+,*,a);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf(" "))
printf("%2d",*elem(a,i,j));
printf("\n");
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf(" "))
printf("%2d",*elem(b,i,j));
printf("\n");

free(a);
free(b);
}

return 0;
}

[asetof...@gmail.com]

luser droog wrote:"
perform a (2D) matrix multiplication upon x and y
using operations F and G.
Z[i,j] = F/ X[i,*] G Y'[j,*]
*/
#define matmul(X,F,G,Y) (matmul)(X,*#F,*#G,Y)
arr (matmul)(arr x, char f, char g, arr y){
int i,j;
arr z=array(x->dims[0],y->dims[1]);

y=clone(y);
transpose2(y);
for (i=0; i<x->dims[0]; i++){
for (j=0; j<y->dims[0]; j++){
arr xs = slice(x,i);
arr ys = slice(y,j);
*elem(z,i,j)=(reduce)(f,(binop)(xs,g,ys));
free(xs);
free(ys);
}
}

free(y);
return z;
}"

A mxn nxk multiplication,
can be defined without the use
of new matrix to free o to traspose
etc
In C++ Something as
i32 mul(matrix<T>& a, matrix<T>& b)
3 matrices a, b and *this as result
Gaming with their index
the loop result something as
for(i=0; i<i0; ++i)
{for(j=0; j<i1; ++j)
for(k=0,s=0u; k<?; ++k)
{s=s+a.e(i,k)*b.e(k,j);}
*this.e(i,j)=s;
}

[asetof...@gmail.com]

[Ike Naar]

On 2015-05-11, luser droog <luser...@gmail.com> wrote:
> void transpose2(arr a){
> int t;
> //if (a->rank != 2) error();
> t = a->dims[0]; a->dims[0] = a->dims[1]; a->dims[1] = t;
> t = a->weight[0]; a->weight[0] = a->weight[1]; a->weight[1] = t;

I'm not sure if the weight array is handled correctly.
Suppose arr a has rank=2, dims={3,4}, weight={4,1}.
Then you want the transposed array to have rank=2, dims={4,3}, weight={3,1},
but not rank=2, dims={4,3}, weight={1,4}.
In other words, the elements of the transposed weight array cannot be obtained
by simply swapping the elements of the original weight array.

> }

[Ike Naar]

Please disregard this post.
The code is okay, I mistook a neat trick for a mistake.
Sorry about that.

[asetof...@gmail.com]

But use of complex macros
with argument can be not good
the time one has to debug...

The max complexity in macros
I have to debug
it is the macro
#define G(a, b) if(a)goto b
(the only macro with argument
I remember I use)
that not show problems in debug.
I have to say the compiler traslate
that in an easy neat way

[luser droog]

At the risk of having some of my own code
shoved in my face, I daresay that G macro
makes things just about impossible to read.
Perhaps with the aid of debugger, but, ...

I have tried to read your code posts to
try to get more ideas, but it's only a
slight exaggeration to say that all I
could understand was "..." (which I
promptly borrowed, so thank you).

I think your G macro is precisely the sort
of jumping that Dijkstra was talking about.
Yes, it's tremendously powerful! Yes, it
can simulate any other control structure
you could possibly want, including the
basic suite. But ... It makes everything
look the same. It's easier to visually
distinguish "if" from "while" than
"G(??)...G(???)...G(??)" from "G(?)...G(?)".

[luser droog]

belated pithy punchline:

That's what the compiler is for!
If you want machine code you know where to find it.

--
sure, if they put the melody back in...

[asetof...@gmail.com]

Using indentation I not see all the
same

For example loop:
l0: Ist1; istr2; G(cond,l0);

Goto used for condition error
G(condErr, err); G(condErr1, err);


Goto used as if()
G(!cond, l0); istr; istr;
l0:

[asetof...@gmail.com]

#include <stdio.h>
#include <stdlib.h>

#include <stdarg.h>
#include <string.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

u32 PBadAlloc; /* if some allocation error for all permutation here 1 */

template<class T> class permutazioni{
public:

permutazioni(u32 a)
{u32 i;
n=0; arr=0; c=0; cnt=0; ji=0;
if(a==0||a>0xFF) {PBadAlloc=1; R;}
arr=(T*)malloc( a*sizeof(T) );
if(arr==0){PBadAlloc=1; R;}
c=(u32*)malloc( a*sizeof(u32) );
if(c==0){free(arr); PBadAlloc=1; arr=0; R;}
n=a;
F(i=0; i<n; ++i){c[i]=0; arr[i]=i;}
}

~permutazioni(){free(c); free(arr);}

u32 swp(T* a, u32 i, u32 j)
{T x;
if(i>=n||j>=n) R -1;
if(i==j) R 0;
x=a[i]; a[i]=a[j]; a[j]=x; ++cnt;
R 0;
}

/* from Sedgewick, result in arr [3!=6] */
u32 nextP(u32 a)
{u32 k;

if(a==0)
{F(k=0; k<n; ++k) c[k]=0;
ji=0; cnt=0;
R 1;
}
F( ; ji<n; )
{if(c[ji]<ji)
{swp(arr, ((ji+1)%2? 0: c[ji]) , ji);
++c[ji]; ji=1; // can be ok ji=0; too
R 1; // ok
}
else {c[ji++]=0; }
}
R 0; // for end
}

void print(void)
{u32 i;
F(i=0; i<n; ++i)
if(i==0 ) cout<<"[ "<< arr[i]<<" ";
else if(i==n-1) cout<<arr[i]<<" ] ";
else cout<<arr[i]<<" ";
}

i32 e(void){R PBadAlloc;}

T& operator[](i32 i){R arr[i];}

u32 n; /* number of element index 0..n-1 */
T *arr; /* where they are */
u32 *c; /* index array for find next permutation*/
u32 ji; /* index for next permutation: nextP() */
u32 cnt; /* cnt%2?-1:1 sign of permutation nextP() return*/
};

// functions for detect u32 or i32 overflow
u32 omul(u32 a, u32 b)
{asm{mov eax, a
mov edx, b
mul edx
cmp edx, 0
je lab
mov eax, -1
lab:
}
R _EAX;
}

u32 oadd(u32 a, u32 b)
{asm{mov eax, a
mov edx, b
add eax, edx
jnc lab
mov eax, -1
lab:
}
R _EAX;
}

u32 osub(u32 a, u32 b)
{asm{mov eax, a
mov edx, b
sub eax, edx
jnc lab
mov eax, -1
lab:
}
R _EAX;
}

i32 omul(i32 a, i32 b)
{asm{mov eax, a
mov edx, b
imul edx
jno lab
mov eax, 0x80000000
lab:
}
R _EAX;
}


i32 oadd(i32 a, i32 b)
{asm{mov eax, a
mov edx, b
add eax, edx
jno lab
mov eax, 0x80000000
lab:
}
R _EAX;
}

i32 osub(i32 a, i32 b)
{asm{mov eax, a
mov edx, b
sub eax, edx
jno lab
mov eax, 0x80000000
lab:
}
R _EAX;
}

// 32 bit unsigned and signed with overflow
class ou32{
public:
ou32(){v=0;}
ou32 operator=( ou32 a){v=a.v; R *this;}
ou32 operator*=(ou32 a){v=omul(v, a.v); R *this;}
ou32 operator+=(ou32 a){v=oadd(v, a.v); R *this;}
ou32 operator-=(ou32 a){v=osub(v, a.v); R *this;}
ou32 operator/=(ou32 a){if(v==-1||a.v==-1) v=-1;
else v=v/a.v;
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, ou32 b){ou32 r; r.v=oadd(a.v, b.v); R r;}
friend ou32 operator-(ou32 a, ou32 b){ou32 r; r.v=osub(a.v, b.v); R r;}
friend ou32 operator/(ou32 a, ou32 b)
{ou32 r;
if(a.v==-1||b.v==-1) r.v=-1;
else r.v=a.v/b.v;
R r;
}
friend ostream& operator<<(ostream& ostr, ou32 a)
{if( a.v==-1 ) ostr<<"ErrElm"; else ostr<<a.v; R ostr;}
friend i32 operator==(ou32 a, ou32 b){R a.v==b.v;}
friend i32 operator!=(ou32 a, ou32 b){R a.v==b.v;}
friend i32 operator<=(ou32 a, ou32 b){R a.v<=b.v;}
friend i32 operator>=(ou32 a, ou32 b){R a.v>=b.v;}
friend i32 operator<( ou32 a, ou32 b){R a.v< b.v;}
friend i32 operator>( ou32 a, ou32 b){R a.v> b.v;}

u32 v;
};

class oi32{
public:
oi32(){v=0;}
operator i32() const{R v;}
oi32 operator=(i32 h){v=h; R *this;}
oi32 operator=(oi32 a){v=a.v; R *this;}
oi32 operator*=(oi32 a){v=omul(v, a.v); R *this;}
oi32 operator+=(oi32 a){v=oadd(v, a.v); R *this;}
oi32 operator-=(oi32 a){v=osub(v, a.v); R *this;}
oi32 operator/=(oi32 a)
{if(v==0x80000000||a.v==0x80000000) v=0x80000000;
else v=v/a.v;
R *this;
}
oi32 operator++(){v=(v==0x80000000? 0x80000000: ++v); R *this;}
oi32 operator++(int){oi32 t=*this; if(v!=0x80000000) ++v; R t;}
oi32 operator--(){v=(v==0x80000000? 0x80000000: --v); R *this;}
oi32 operator--(int){oi32 t=*this; if(v!=0x80000000) --v; R t;}
friend oi32 operator*(oi32 a, oi32 b){oi32 r; r.v=omul(a.v, b.v); R r;}
friend oi32 operator+(oi32 a, oi32 b){oi32 r; r.v=oadd(a.v, b.v); R r;}
friend oi32 operator-(oi32 a, oi32 b){oi32 r; r.v=osub(a.v, b.v); R r;}
friend oi32 operator/(oi32 a, oi32 b)
{oi32 r;
if(a.v==0x80000000||b.v==0x80000000) r.v=0x80000000;
else r.v=a.v/b.v;
R r;
}
friend ostream& operator<<(ostream& ostr, oi32 a)
{if( a.v==0x80000000 ) ostr<<"ErrElm"; else ostr<<a.v; R ostr;}
friend i32 operator==(oi32 a, oi32 b){R a.v==b.v;}
friend i32 operator!=(oi32 a, oi32 b){R a.v==b.v;}
friend i32 operator<=(oi32 a, oi32 b){R a.v<=b.v;}
friend i32 operator>=(oi32 a, oi32 b){R a.v>=b.v;}
friend i32 operator<( oi32 a, oi32 b){R a.v< b.v;}
friend i32 operator>( oi32 a, oi32 b){R a.v> b.v;}

i32 v;
};

u32 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: MtrBadAlloc=1; R -1; // for error

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
i32 ov(i32 a){R 0;}
i32 ov(d64 a){R a==NAN||a==+INF||a==-INF ? 1: 0; } // overflow con NAN
i32 ov(ou32 a){R a.v==-1? 1: 0; } // overflow con -1 range: 0..0FFFFFFFE
i32 ov(oi32 a){R a.v==0x80000000? 1: 0; } // overflow con 0x80000000

Matrix<T>& operator=(Matrix<T>& a)
{u32 i;
if(a.mtr== 0)
{
l0: *Led(mtr)=1; G(1,z); /* errore trovato e riportato */
}

G(*Led(a.mtr), l0);

if(a.mtr!=mtr)
{// deve essere delle stesse dimensioni e indici massimi
G(a.mtr[0]!=mtr[0], l0);
F(i=1;i<=mtr[0];++i)
G(a.mtr[i]!=mtr[i], l0);
G(*Led(a.mtr)!=0, l0); G(*SZe(mtr)!=*SZe(a.mtr), 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
G(a.mtr[0]!=mtr[0], l0);
F(i=1;i<=mtr[0];++i)
G(a.mtr[i]!=mtr[i], l0);
G(*Led(a.mtr)!=0, 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
G(a.mtr[0]!=mtr[0], l0);
F(i=1;i<=mtr[0];++i)
G(a.mtr[i]!=mtr[i], l0);
G(*Led(a.mtr)!=0, 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*=(T& a)

{u32 i;
T *p1;
u8 *pmax;

G(*Led(mtr), z);

p1=(T*) Data(mtr);
i=*SZ(mtr); pmax=(u8*) p1; pmax+=i;
F(; p1<pmax ;++p1)
{*p1 = *p1 * a;

if(ov(*p1)) {*Led(mtr)=1; G(1,z);}

}
z:
R *this;
}

Matrix<T>& operator*=(T a)
{u32 i;
T *p1;
u8 *pmax;

G(*Led(mtr), z);

p1=(T*) Data(mtr);
i=*SZ(mtr); pmax=(u8*) p1; pmax+=i;
F(; p1<pmax ;++p1)
{*p1 = *p1 * a;

if(ov(*p1)) {*Led(mtr)=1; G(1,z);}
}
z:
R *this;
}

// *this=a*b
// ove *this e' come a mxn e b e' nxj
i32 mul(Matrix<T>& a, Matrix<T>& b)
{u32 i, j, z, i0, i1;
T a1, a2, s, t1, t2;

G( *Led(mtr) , l1);
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; }
if( *SZe(mtr) != *SZe(a.mtr) )
G(1, l0);
if(a[1]!=b[0]||(*this)[0]!=a[0]||(*this)[1]!=a[1])
G(1, l0);
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;
}
}
R 0;
}

i32 sub(Matrix<T>& a, Matrix<T>& b)
{u32 i;
T *p1,*p2, p3;
u8 *pmax;

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
G(a.mtr[0]!=mtr[0], l0);

G(b.mtr[0]!=mtr[0], l0);

F(i=1;i<=mtr[0];++i)

G(a.mtr[i]!=mtr[i]||b.mtr[i]!=mtr[i], l0);
G(*Led(a.mtr)!=0, l0);
G(*Led(b.mtr)!=0, 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 add(Matrix<T>& a, Matrix<T>& b)
{u32 i;
T *p1,*p2, p3;
u8 *pmax;

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);

G(a.mtr[0]!=mtr[0], l0);

G(b.mtr[0]!=mtr[0], l0);

F(i=1;i<=mtr[0];++i)

G(a.mtr[i]!=mtr[i]||b.mtr[i]!=mtr[i], l0);
G(*Led(a.mtr)!=0, l0);
G(*Led(b.mtr)!=0, 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;
if(a.mtr==0)
{
l0: *Led(mtr)=1; l1: R -1;
}
G(*Led( mtr), l1);
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[3]=a.mtr[4]; mtr[4]=a.mtr[3];
R 0;
}

i32 transpose(void)
{u32 t;
if(mtr[0]!=2)
{
l0: *Led(mtr)=1; l1: R -1;
}
G( *Led(mtr), l1 );
t=mtr[1]; mtr[1]=mtr[2]; mtr[2]=t;
t=mtr[3]; mtr[3]=mtr[4]; mtr[4]=t;
R 0;
}

T det(void)
{u32 t, i1, k;
permutazioni<u32> pe(mtr[1]);
T x, p, s;
s=-1;
if(mtr[0]!=2u)
{
l0: *Led(mtr)=1;
l1: R s;
}
G( *Led(mtr), l1);
G(mtr[1]!=mtr[2], l0);
i1=mtr[1]; s=0;
la: F(k=0, p=1; k<i1; ++k)
p=p*(*this).e(k, pe[k]);
G(ov(p),l0);
if(pe.cnt%2!=0) {s=s-p; // cout<<"-"<<(i32)p;
}
else {s=s+p; // cout<<"+"<<(i32)p;
}
G(ov(s),l0);
G( (t=pe.nextP(1))==-1, l0);
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;
p1=(T*) Data(mtr); p2=a;

i=*SZ(mtr); pmax=(u8*) p1; pmax+=i;
F(; p1<pmax ;++p1, ++p2)

*p1=*p2;

}

/* ritorna il sup dell'indice consentito */
/* indici = 0, 1, 2 ,3 ,4 */

i32 operator[](i32 i)

{if(mtr==0||mtr[0]<=i) R 0;
R mtr[i+1];
}
int print(void){PrntMtrx(mtr);}

u32 *mtr;
};

u32 f(void)
{Matrix<u32> m1(2,2),m2(2,3),m6(2,2),m7(2,2),m8(1,2),m9(1,2);
Matrix<d64> m3(2,2);
Matrix<ou32> m10(2,2);
Matrix<u8> m4(1,7), m5(7,1);
Matrix<oi32> m12(4,4), m13(3,3);
u32 d, i, j;
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;

d=m1.det();

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");
P("M5:\n"); m5.print(); P("\n\n");

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);
m7=m1;
P("M7: m7=m1\n"); m7.print(); P("\n\n");

m1*=10;
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);
R 0;
}

i32 g(void)
{u32 i, r, j;
permutazioni<u32> m(3);
permutazioni<u8> v(4);
char *ptr="mare";

if(m.e()) R 0;
F(i=0; i<m.n; ++i)
m[i]=i;
F(;;){m.print();
G( (r=m.nextP(1))==-1 || r==0, m1);
}
m1:
F(i=0; i<v.n; ++i)
v[i]=ptr[i];

F(;;){v.print();
G( (r=v.nextP(1))==-1 || r==0, la);
}
la:
P("\n");
R 0;
}

int main(void)
{if( sizeof(u32)!=4 || sizeof(u16)!=2 )
{P("Basic type error\n"); R 0;}
f();
R 0;
}

------------
M1:
1 2
3 4


m1.det()=-2

M2:
0 1 2
10 11 12


M3:
0.00 1.00
10.23 11.23


M4:
0 1 2 3 4 5 6


M5:
0
10
20
30
40
50
60

M6:
5 6
7 8


M7=m1*m6:
19 22
43 50


M8
1 2


M9=m8*m1
7 10


M12
2 9 9 4
2 4294967293 12 8
4 8 3 4294967291
1 2 6 4


DET(M12)=147
M7: m7=m1
1 2
3 4


M1=M1*10:
10 20
30 40


M1+=M1:
20 40
60 80


M13:

1 2 3
4 5 6
7 8 9

DET(M13)=0

[asetof...@gmail.com]

I wrote:"
T det(void)
{u32 t, i1, k;
permutazioni<u32> pe(mtr[1]);
T x, p, s;
s=-1;
if(mtr[0]!=2u)
{
l0: *Led(mtr)=1;
l1: R s;
}
G( *Led(mtr), l1);
G(mtr[1]!=mtr[2], l0);
i1=mtr[1]; s=0;
la: F(k=0, p=1; k<i1; ++k)
p=p*(*this).e(k, pe[k]);
G(ov(p),l0);
if(pe.cnt%2!=0) {s=s-p; // cout<<"-"<<(i32)p;
}
else {s=s+p; // cout<<"+"<<(i32)p;
}
G(ov(s),l0);
G( (t=pe.nextP(1))==-1, l0);
G( t==1, la);
R s;
}"

This above it is bugged
Perhaps this is better:

T det(void)
{u32 t, i1, k;
permutazioni<u32> pe(mtr[1]);

T rtn, p, s;
rtn=returnOverflowTelement(rtn);
if(mtr[0]!=2u)
{
l0: R rtn;
}
G( *Led(mtr), l0);

[luser droog]

On Monday, May 11, 2015 at 9:01:47 PM UTC-5, asetof...@gmail.com wrote:
> I wrote:"
> A mxn nxk multiplication,
> can be defined without the use
> of new matrix to free o to traspose
> etc
> In C++ Something as
> i32 mul(matrix<T>& a, matrix<T>& b)
> 3 matrices a, b and *this as result
> Gaming with their index
> the loop result something as

> for(i=0; i<i0; ++i)
> {for(j=0; j<i1; ++j)
> {for(k=0,s=0u; k<?; ++k)
> {s=s+a.e(i,k)*b.e(k,j);}
> *this.e(i,j)=s;
> }
> }

Took a little more time and energy than anticipated,
but finally here's the cake. matmul() now implements
the generalized Inner Product for arguments of
arbitrary, compatible dimensions.

Getting a little long for this format. Probably the
last post I'll make here. This file is also available
at https://github.com/luser-dr00g/inca/blob/master/arrind.c

//arrind.c - Arrays, indirect.

#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "ppnarg.h" //https://github.com/luser-dr00g/inca/ppnarg.h

typedef struct arr {
int rank;
int *dims;
int *weight;
int *data;
} *arr;

/* multiply together rank integers in dims array
*/
int productdims(int rank, int *dims){
int i,z=1;
for(i=0; i<rank; i++)
z *= dims[i];
return z;
}

arr arraya(int rank, int *dims){

int datasz;
int i;
int x;
arr z;

datasz=productdims(rank,dims);
z=malloc(sizeof(struct arr)
+ (rank+rank+datasz)*sizeof(int));

z->rank = rank;
z->dims = (int*)(((char*)z) + sizeof(struct arr));
z->weight = z->dims + rank;
z->data = z->weight + rank;
memmove(z->dims,dims,rank*sizeof(int));
for(x=1, i=rank-1; i>=0; i--){
z->weight[i] = x;
x *= z->dims[i];
}

return z;
}

void loaddimsv(int rank, int dims[], va_list ap){
int i;

for (i=0; i<rank; i++){
dims[i]=va_arg(ap,int);
}
}

/* create a new array with specified dimensions */

#define array(...) (array)(PP_NARG(__VA_ARGS__),__VA_ARGS__)

/* create a new array with specified rank and dimensions
*/

arr (array)(int rank, ...){

va_list ap;
//int *dims=calloc(rank,sizeof(int));
int dims[rank];

int i;
int x;
arr z;

va_start(ap,rank);

loaddimsv(rank,dims,ap);
va_end(ap);

z = arraya(rank,dims);
//free(dims);
return z;
}

/* create an array header to access existing data in multidimensional layout */
arr cast(int *data, int rank, ...){
va_list ap;
int dims[rank];

int i;
int x;
arr z;

va_start(ap,rank);

loaddimsv(rank,dims,ap);
va_end(ap);

z=malloc(sizeof(struct arr)
+ (rank+rank)*sizeof(int));

z->rank = rank;

z->dims = (int *)(((char *)z) + sizeof(struct arr));

z->weight = z->dims + rank;

z->data = data;

memmove(z->dims,dims,rank*sizeof(int));
for(x=1, i=rank-1; i>=0; i--){
z->weight[i] = x;
x *= z->dims[i];
}

return z;
}

/* create a new array which shares the data of an existing array
*/
arr clone(arr a){
arr z=malloc(sizeof(struct arr)
+ (a->rank+a->rank)*sizeof(int));
z->rank = a->rank;
z->dims = (int*)(((char*)z) + sizeof(struct arr));
z->weight = z->dims + z->rank;
z->data = a->data;
memmove(z->dims,a->dims,z->rank*sizeof(int));
memmove(z->weight,a->weight,z->rank*sizeof(int));
return z;
}

/* exchange the leftmost two dimensions (only two in 2D)
*/
void transpose2(arr a){
int t;
//if (a->rank != 2) error();
t = a->dims[0]; a->dims[0] = a->dims[1]; a->dims[1] = t;
t = a->weight[0]; a->weight[0] = a->weight[1]; a->weight[1] = t;
}

/* rotate dims and weights according to sign and magnitude of shift
transpose(1,a)==transpose(-1,a)==transpose2(a) for 2D
*/
void transpose(int shift, arr a){
int i;
int t;
while(shift){
if (shift>0){
t=a->dims[0];
for (i=1; i<a->rank; i++)
a->dims[i-1]=a->dims[i];
a->dims[a->rank-1]=t;
t=a->weight[0];
for (i=1; i<a->rank; i++)
a->weight[i-1]=a->weight[i];
a->weight[a->rank-1]=t;
--shift;
} else {
t=a->dims[a->rank-1];
for (i=a->rank-2; i>=0; i--)
a->dims[i+1]=a->dims[i];
a->dims[0]=t;
t=a->weight[a->rank-1];
for (i=a->rank-2; i>=0; i--)
a->weight[i+1]=a->weight[i];
a->weight[0]=t;
++shift;

/* compute vector index list for ravel index ind */
int *vector_index(int ind, int *dims, int n, int *vec){
int i,t=ind, *z=vec;

for (i=0; i<n; i++){

z[n-1-i] = t % dims[n-1-i];
t /= dims[n-1-i];
}
return z;
}

/* compute ravel index for vector index list */
int ravel_index(int *vec, int *dims, int n){
int i,z=*vec;
for (i=0; i<n-1;i++){
z*=dims[i+1];
z+=vec[i+1];
}
return z;
}

/* create a vector of all elements of x followed by all elements of y */
arr catv(arr x, arr y){
int xsz = productdims(x->rank,x->dims);
int ysz = productdims(y->rank,y->dims);
int datasz = xsz + ysz;
arr z=array(datasz);
int scratch[x->rank+y->rank];
int i;
for (i=0; i<xsz; i++)
*elem(z,i) = *elema(x,vector_index(i,x->dims,x->rank,scratch));
for (i=0; i<ysz; i++)
*elem(z,xsz+i) = *elema(y,vector_index(i,y->dims,y->rank,scratch));

return z;
}

/* create a (contiguous) copy of a (not necessarily contiguous)
existing array
*/
arr copy(arr a){
int datasz = productdims(a->rank,a->dims);
arr z=malloc(sizeof(struct arr)
+ (a->rank+a->rank+datasz)*sizeof(int));
int i;
int x;
int ind[a->rank];

z->rank = a->rank;
z->dims = (int*)(((char*)z) + sizeof(struct arr));
z->weight = z->dims + z->rank;
z->data = z->weight + z->rank;
memmove(z->dims,a->dims,z->rank*sizeof(int));
for(x=1, i=z->rank-1; i>=0; i--){
z->weight[i] = x;
x *= z->dims[i];
}

for (i=0;i<datasz;i++){
int j;

#if 0

int idx = i;
for (j=z->rank-1; j>=0; j--) { //generate index
ind[j] = idx % z->dims[j];
idx /= z->dims[j];
}

#endif
vector_index(i,z->dims,z->rank,ind);

/* perform a (2D) matrix multiplication upon rows of x and columns of y

using operations F and G.

Z = X F.G Y

Z[i,j] = F/ X[i,*] G Y'[j,*]

more generally,
perform an inner product on arguments of compatible dimension.
Z = X[A;B;C;D;E;F] +.* Y[G;H;I;J;K] |(F = G)
Z[A;B;C;D;E;H;I;J;K] = +/ X[A;B;C;D;E;*] * Y[*;H;I;J;K]

*/
#define matmul(X,F,G,Y) (matmul)(X,*#F,*#G,Y)
arr (matmul)(arr x, char f, char g, arr y){
int i,j;

//arr z=array(x->dims[0],y->dims[y->rank-1]);
arr xdims = cast(x->dims,1,x->rank);
arr ydims = cast(y->dims,1,y->rank);
xdims->dims[0]--;
ydims->dims[0]--;
ydims->data++;
arr z=arraya(x->rank+y->rank-2,catv(xdims,ydims)->data);
int datasz = productdims(z->rank,z->dims);
int k=y->dims[0];
arr xs = array(k);
arr ys = array(k);

y=clone(y);
transpose(1,y);

for (i=0; i<datasz; i++){

int idx[z->rank];
vector_index(i,z->dims,z->rank,idx);
int *xdex=idx;
int *ydex=idx+x->rank;
memmove(ydex,ydex-1,y->rank);
for (j=0; j<k; j++){
xdex[x->rank-1]=j;
xs->data[j] = *elema(x,xdex);
}
for (j=0; j<k; j++){
ydex[y->rank-1]=j;
ys->data[j] = *elema(y,ydex);
}
z->data[i] = (reduce)(f,(binop)(xs,g,ys));
}

free(y);
return z;
}

arr iota(int n){
arr z = array(n);

int i;
for (i=0;i<n;i++)

*elem(z,i)=i;
return z;
}

int main(){

#ifdef TEST_BASIC

#ifdef TEST_MATMUL

for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf(" "))

printf("%3d",*elem(a,i,j));
printf("\n");

b = matmul(a,+,*,a);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf(" "))

printf("%3d",*elem(b,i,j));

printf("\n");

free(a);
free(b);
}

{
int i,j,k,l,n=2;
arr a=iota(n*n*n);
arr b=cast(a->data,3,n,n,n);
arr c=matmul(b,+,*,b);
printf("%d\n",c->rank);
for (i=0; i<c->rank; i++,printf(" "))
printf("%d",c->dims[i]);
printf("\n");
for (i=0; i<c->dims[0]; i++,printf("\n"))
for (j=0; j<c->dims[1]; j++,printf("\n"))
for (k=0; k<c->dims[2]; k++,printf("\n"))
for (l=0; l<c->dims[3]; l++,printf(" "))
printf("%3d", *elem(c,i,j,k,l));

printf("\n");
free(a);
free(b);
free(c);
}

#endif

#ifdef TEST_CAST
{ /* testing cast() */
int i,j,k,n=3;
int q[n][n][n];
arr a=cast((int *)q,3,n,n,n);

for (i=0;i<n;i++)
for (j=0;j<n;j++)
for (k=0;k<n;k++)
*elem(a,i,j,k) = n*n*n - ((i*n+j)*n+k);
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf("\n"))
for (k=0;k<n;k++,printf(" "))
printf("%2d",*elem(a,i,j,k));
printf("\n");
for (i=0;i<n;i++,printf("\n"))
for (j=0;j<n;j++,printf("\n"))
for (k=0;k<n;k++,printf(" "))

printf("%2d", q[i][j][k]);
printf("\n");

free(a);
}
#endif

#ifdef TEST_IOTA
{
arr a=iota(12);
arr b=catv(a,a);
int i;
for (i=0;i<b->dims[0];i++,printf(" "))
printf("%2d", *elem(b,i));

arr c=cast(b->data,3,2,3,4);
int j,k;
for (i=0; i<2; i++,printf("\n"))
for (j=0; j<3; j++,printf("\n"))
for (k=0; k<4; k++,printf(" "))
printf("%2d", *elem(c,i,j,k));

free(a);
free(b);
free(c);
}
#endif

return 0;
}

--
eof

[asetof...@gmail.com]

here the traspose of
a=0 1 2 3 4 5 6
would be b=

0
1
2
3
4
5
6

for a 'rank'=2 i0=1 i1=7 w0=7 w1=1
for b 'rank'=2 i0=7 i1=1 w0=1 w1=1

So b here can not be obtained from
a just swap i0<->i1 and w0<->w1

[asetof...@gmail.com]

Here this
1 2 3 4
5 6 7 8
has rank=2 i0=2 i1=4 w0=4 w1=1

5 1
6 2
7 3
8 4
has rank=2 i0=4 i1=2 w0=2 w1=1
And would be one the transpose
of the other...
I think for traspose swap i0 i1
and from them recalculate w0 and w1.
This would be not sufficient because
for me one has to swap elements in
array too...

[asetof...@gmail.com]

Here this
1 2 3 4
5 6 7 8
has rank=2 i0=2 i1=4 w0=4 w1=1

1 5
2 6
3 7
4 8

has rank=2 i0=4 i1=2 w0=2 w1=1
And would be one the transpose
of the other...
I think for traspose swap i0 i1
and from them recalculate w0 and w1.
This would be not sufficient because
for me one has to swap elements in
array too...

The linear elements of 1 are

1 2 3 4 5 6 7 8

The linear elements of 2 are
1 5 2 6 3 7 4 8

So in this case for to have the
right matrix I would change
linear array (or data) for the
transpose matrix too
not only i0 etc

[luser droog]

I think it still works for this case. Remember that
this is not a normal 'transpose' like everybody else
does. We don't rearrange the data at all. This is more
properly called a 'dope-vector pseudo-transpose'.
Instead of changing the data, we just change the
constants in the access formula, rearranging the
coefficients of the polynomial. In a real sense, this
is just an application of the commutativity and
associativity of addition.

So for concreteness, assume the data is arranged
sequentially starting at hypothetical address 500.

500: 0
501: 1
502: 2
503: 3
504: 4
505: 5
506: 6

if a is rank 2, dims {1, 7), weight (7, 1), then the
sum of the indices multiplied by the associated weights
added to the initial pointer (500) yield the appropriate
addresses for each element
a[0][0] == *(500+0*7+0*1)
a[0][1] == *(500+0*7+1*1)
a[0][2] == *(500+0*7+2*1)
a[0][3] == *(500+0*7+3*1)
a[0][4] == *(500+0*7+4*1)
a[0][5] == *(500+0*7+5*1)
a[0][6] == *(500+0*7+6*1)

So the dope-vector pseudo-transpose rearranges the
weights and dimensions to match the new ordering of
indices, but the sum remains the same. The initial
pointer remains the same. The data does not move.

b[0][0] == *(500+0*1+0*7)
b[1][0] == *(500+1*1+0*7)
b[2][0] == *(500+2*1+0*7)
b[3][0] == *(500+3*1+0*7)
b[4][0] == *(500+4*1+0*7)
b[5][0] == *(500+5*1+0*7)
b[6][0] == *(500+6*1+0*7)

[luser droog]

On Tuesday, May 19, 2015 at 1:16:30 PM UTC-5, asetof...@gmail.com wrote:
> Here this
> 1 2 3 4
> 5 6 7 8
> has rank=2 i0=2 i1=4 w0=4 w1=1
>
> 5 1
> 6 2
> 7 3
> 8 4
> has rank=2 i0=4 i1=2 w0=2 w1=1
> And would be one the transpose
> of the other...

Here you've done a different kind of transpose, about a
different axis. Or you could consider it a transpose +
reversal. A "simple" transpose would yield

1 5
2 6
3 7
4 8

> I think for traspose swap i0 i1
> and from them recalculate w0 and w1.
> This would be not sufficient because
> for me one has to swap elements in
> array too...

I think this kind of transpose can also be accommodated
by using a negative weight, but I haven't fully thought
it through. You might also need to adjust the pointer
to account for the "new" first element.

[asetof...@gmail.com]

Transpose here is ok...
It seems code write itself...

800 lines in C++ like
Matrix on numeric type
operator defined
=,+,-,*, determinant, transpose
print

base type build with overflow
that matrix class can be used
detecting overflow
To do: Inverse, <<, >>, ==
'complementi algebrici'

[asetof...@gmail.com]

luser droog wrote:"
I think it still works for this case. Remember that
this is not a normal 'transpose' like everybody else
does. We don't rearrange the data at all. This is more
properly called a 'dope-vector pseudo-transpose'.
Instead of changing the data, we just change the
constants in the access formula, rearranging the
coefficients of the polynomial. In a real sense, this
is just an application of the commutativity and
associativity of addition."

I rearrange data instead on traspose
matrices.
Yes could be computationally not
good but
I'm sure in this way that i j index
near has a near collocation in the
memory matrix linear layout
And all matrices as
A11 A12 ... A1n
...............
Am1 Am2.....Amn
follow this linear memory configuration
A11,A12,..A1n,A2n....Am1,Am2,...Amn

[asetof...@gmail.com]

I wrote:"
I rearrange data instead on traspose
matrices.
Yes could be computationally not
good but
I'm sure in this way that i j index
near has a near collocation in the
memory matrix linear layout
And all matrices as
A11 A12 ... A1n
...............
Am1 Am2.....Amn
follow this linear memory configuration
A11,A12,..A1n,A2n....Am1,Am2,...Amn
"

Yes I can make some errors on that

[asetof...@gmail.com]

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;

[asetof...@gmail.com]

I wrote:"
adet=2; adet/=2;
G(adet==0, l0); /* il tipo non ammette inversa */"

This would be:
adet=1; adet/=2;
G(adet==0, l0);
It would be some float Point type etc

[luser droog]

On Wednesday, May 6, 2015 at 1:18:25 AM UTC-5, luser droog wrote:
> I've been thinking further about the arbitrary-dimensional array
> structure (and am still stunned by the lackluster response so far),
> and I think a suitable coup-de-grace for demonstrating the value
> of indirection is to tease-out yet another data member for this
> structure. This ties everything back to the 1962 APL book, the
> chapter on Representation, which has kept me up so many, many nights.
>
> A dynamic array object is defined by its rank, list of dimension sizes,
> a pointer to the start of the data, AND a weighting vector which for
> simple uses merely caches the multiplications from the array-access
> formula but it affords us much greater flexibility as (I hope) we
> shall see...
>

Posted a new write-up on SO
http://stackoverflow.com/questions/30409991/use-a-dope-vector-to-access-arbitrary-slices-of-a-multidimensional-array

--
y'all know how to upvote, right?

[asetof...@gmail.com]

And how do you print the values
(the elements )
of one 3 index Matrix?
or 4 index matrix?

I had seen such matrices only in
tensors mathematics
That made sense after some index is lost
thru multiplication with a
vector or some other matrix

[asetof...@gmail.com]

So where are useful these >2 dimensions matrices?

[luser droog]

Alright, I came up with some code to print them.
Until now, I've just been printing them "manually"
while testing. And the printing functions in all
the versions of inca are pretty terrible and limited.

I tried to write this recursively/inductively, so it
*ought* to be correct, without extensive testing,
because it *looks* correct. I think. If I did it right.

And you're right about the relation to tensor
mathematics. Iverson borrowed lots of ideas there.

other functions at
https://github.com/luser-dr00g/inca/blob/master/arrind.c

int iscontiguous(arr a){
int i,x;
for (i=a->rank-1,x=1; i>=0; i--){
if (a->weight[i] != x)
return 0;
x *= a->dims[i];
}
return 1;
}

void print(arr a, int width){
int i;
int maxwidth;
int freecopy = 0;

if (width){
maxwidth=width;
} else {

int datasz = productdims(a->rank,a->dims);

if (!iscontiguous(a)) {
a = copy(a);
freecopy = 1;
}

maxwidth=0;

for (i=0; i<datasz; i++){

int size = snprintf(NULL,0,"%d",a->data[i]);
if (size > maxwidth)
maxwidth = size;
}
}

switch(a->rank){
case 0: printf("%*d\n", maxwidth, *a->data);
break;
case 1: for (i=0; i<a->dims[0]; i++)
printf("%*d ", maxwidth, *elem(a,i));
break;
default:
for (i=0; i<a->dims[0]; i++){
arr as = slice(a,i);
print(as,maxwidth);
printf("\n");
free(as);
}
break;
}

if (freecopy)
free(a);
}


int main() {

#ifdef TEST_PRINT
{
arr a = iota(27);
arr b = casta(a->data, 3, (int[]){3,3,3});
print(b,0);

free(b);
free(a);

a = iota(64);
b = casta(a->data, 4, (int[]){2,2,2,2});
print(b,0);

free(b);
free(a);
}
#endif

return 0;
}


produces:

0 1 2
3 4 5
6 7 8

9 10 11
12 13 14
15 16 17

18 19 20
21 22 23
24 25 26

0 1
2 3

4 5
6 7

8 9
10 11

12 13
14 15

[luser droog]

On Saturday, May 23, 2015 at 11:26:07 PM UTC-5, asetof...@gmail.com wrote:
> So where are useful these >2 dimensions matrices?

I find this to be a very difficult question to answer without some
handwaving. Frankly, wherever you have >1 of these 2d objects that
must be treated as a unit.

But more peosaically,

int seraphim=3,
cherubim=3,
thrones=3,
dominions=3,
powers=3,
authorities=3,
principalities=3,
archangels=3,
angels=3;

arr aleph = array(seraphim,
cherubim,
thrones,
dominions,
powers,
authorities,
principalities,
archangels,
angels);

and eternitie in an houre.

[asetof...@gmail.com]

luser droog wrote:...

Thank you

[asetof...@gmail.com]

// functions for detect u32 or i32 overflow
u32 omul(u32 a, u32 b)
{asm{mov eax, a

cmp eax, -1
je lab
mov edx, b
cmp edx, -1
je m1

mul edx
cmp edx, 0
je lab

m1:

mov eax, -1
lab:
}
R _EAX;
}

u32 oadd(u32 a, u32 b)
{asm{mov eax, a
mov edx, b
add eax, edx
jnc lab
mov eax, -1
lab:
}
R _EAX;
}

u32 osub(u32 a, u32 b)
{asm{mov eax, a

cmp eax, -1
je lab
mov edx, b
cmp edx, -1
je m1

sub eax, edx
jnc lab

m1:

mov eax, -1
lab:
}
R _EAX;
}

i32 omul(i32 a, i32 b)
{asm{mov eax, a

cmp eax, 0x80000000
je lab
mov edx, b
cmp edx, 0x80000000
je m1
imul edx
jno lab
m1:

mov eax, 0x80000000
lab:
}
R _EAX;
}


i32 oadd(i32 a, i32 b)
{asm{mov eax, a

cmp eax, 0x80000000
je lab
mov edx, b
cmp edx, 0x80000000
je m1

add eax, edx
jno lab

m1:

mov eax, 0x80000000
lab:
}
R _EAX;
}

i32 osub(i32 a, i32 b)
{asm{mov eax, a

cmp eax, 0x80000000
je lab
mov edx, b
cmp edx, 0x80000000
je m1

sub eax, edx
jno lab

m1:

[asetof...@gmail.com]

I wrote:"friend ou32 operator*(ou32 a, ou32 b){ou32 r; r.v=omul(a.v, b.v); R r;}"
This should be wrong because here
it return the address in the stack
of r. But as all it can be easily
corrected using one ou32
circular array

[Keith Thompson]

asetof...@gmail.com writes:
> I wrote:"friend ou32 operator*(ou32 a, ou32 b){ou32 r; r.v=omul(a.v, b.v); R r;}"

That's C++, not C.

--
Keith Thompson (The_Other_Keith) ks...@mib.org <http://www.ghoti.net/~kst>
Working, but not speaking, for JetHead Development, Inc.
"We must do something. This is something. Therefore, we must do this."
-- Antony Jay and Jonathan Lynn, "Yes Minister"

[Ian Collins]

Keith Thompson wrote:
> asetof...@gmail.com writes:
>> I wrote:"friend ou32 operator*(ou32 a, ou32 b){ou32 r; r.v=omul(a.v, b.v); R r;}"
>
> That's C++, not C.

I'd call it gobbledegook!

--
Ian Collins

[asetof...@gmail.com]

I wrote:"
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 det(T* result)
Should be the ok as prototype...

How can note G(,) macro it is the
right way for build such functions

[Ian Collins]

asetof...@gmail.com wrote:
> I wrote:"

More gobbledegook?

--
Ian Collins

[mark.b...@gmail.com]

On Monday, 1 June 2015 09:36:14 UTC+1, Ian Collins wrote:
> asetof...@gmail.com wrote:
> > I wrote:"
>
> More gobbledegook?

His approach to C seems to me to be that of someone who really doesn't like the language and wants to punish it.

[luser droog]

On Monday, June 1, 2015 at 3:36:14 AM UTC-5, Ian Collins wrote:
> asetof...@gmail.com wrote:
> > I wrote:"
>
> More gobbledegook?
>

If you like gobbledegook, here's my code all golfed-
up and gook. There's even a little whitespace in there.
A very little.

$ cat arrg.c
#include<stdarg.h>
#include<stdlib.h>
#include<stdio.h>
#include<string.h>
#include"ppnarg.h" //https://github.com/luser-dr00g/inca/blob/master/ppnarg.h
typedef char C; typedef int I;
typedef struct a{I r,*d,*w,*p;}*A;
#define R return
#define DO(n,x){I i=0,N=(n);for(;i<N;++i){x;}}/*loop i=[0..n)*/
I tr(I r,I*d){ I z=1; DO(r,z*=d[i]) R z; } /* productdims */
dw(I r,I*d,I*w){ I x=1; DO(r,w[r-1-i]=x; x*=d[r-1-i]) } /* calculate weights */
A ah(I r,I*d,I sz){ A z=malloc((sizeof*z)+(sizeof(I)*(r+r+sz))); /* array header + sz extra */
z->r=r; z->d=(I*)(((C*)z)+sizeof*z); z->w=z->d+r; z->p=NULL;
memcpy(z->d,d,r*sizeof(I)); dw(r,d,z->w); R z; }
A ara(I r,I*d){ I sz=tr(r,d); A z=ah(r,d,sz); z->p=z->w+r; R z; } /* arraya */
dv(I r,I*d,va_list v){ DO(r,d[i]=va_arg(v,I)) } /* loaddimsv */
A ar(I r,...){ va_list v; I d[r]; va_start(v,r); dv(r,d,v); va_end(v); R ara(r,d); } /* array rm */
#define ar(...)ar(PP_NARG(__VA_ARGS__),__VA_ARGS__)
A caa(I*p,I r,I*d){ A z=ah(r,d,0); z->p=p; R z; } /* casta rm */
A ca(I*p,I r,...){ va_list v; I d[r]; va_start(v,r); dv(r,d,v); va_end(v); R caa((I*)p,r,d); }
#define ca(p,...)ca((I*)p,PP_NARG(__VA_ARGS__),__VA_ARGS__)
A cl(A a){ A z=ah(a->r,a->d,0);
memcpy(z->d,a->d,z->r*sizeof(I));
memcpy(z->w,a->w,z->r*sizeof(I)); z->p=a->p; R z; } /* clone */
I*ela(A a,I*j){ I x=0; DO(a->r, x+=j[i]*a->w[i]) R a->p+x; } /* elema */
I*elv(A a,va_list v){ I j[a->r]; dv(a->r,j,v); R ela(a,j); } /* elemv */
I*el(A a,...){ I*z; va_list v; va_start(v,a); z=elv(a,v); va_end(v); R z; } /* elem */
int *vx(I x,I*d,I r,I*v){ DO(r, v[r-1-i]=x%d[r-1-i]; x/=d[r-1-i]) R v; } /* vector_index */
int rx(I*v,I*d,I r){ I z=*v; DO(r-1, z*=d[i+1]; z+=v[i+1]) R z; } /* ravel_index */
A cp(A a){ I sz=tr(a->r,a->d); A z=ah(a->r,a->d,sz); z->p=z->w+z->r;
I j[z->r]; DO(sz, vx(i,z->d,z->r,j); z->p[i]=*ela(a,j)) R z; } /* copy rm */
#define CASE ;break;case

#define OP(x,y)((((x))*256)+((y)))
#define OPS(_) _("+", '+', 0, +, 0) \
_("*", '*', 0, *, 1) \
_("=", '=', 0, ==, 1) \
_("==",'=','=',==, 1) \
/* f f0 f1 F id */
#define BINOP(f,f0,f1,F,id) CASE OP(f0,f1): *el(z,i) = *el(x,i) F *el(y,i);
A bin(A x,C*f,A y){ A z=cp(x); DO(x->d[0], switch(OP(*f,f[1])){ OPS(BINOP) }) R z; }
#define bin(X,F,Y) bin(X,#F,Y)
#define REDID(f,f0,f1,F,id) CASE OP(f0,f1): z = id;
#define REDOP(f,f0,f1,F,id) CASE OP(f0,f1): z = *el(y,n-2-i) F z;
int red(C*f,A y){ I z; I n=y->d[0];
switch(OP(*f,f[1])){ OPS(REDID) }
if (n){ z=*el(y,n-1);
if (n-1>0) DO(n-1, switch(OP(*f,f[1])){ OPS(REDOP) }) }
R z; }
#define red(F,Y) red(#F,Y)

A xt(A a,I x){ I r=a->r+x; I d[r];
DO(x, d[i]=1) memcpy(d+x,a->d,a->r*sizeof(I)); R caa(a->p,r,d); } /* extend */
A iot(I n){ A z=ar(n); DO(n,*el(z,i)=i) R z; } /* index generator */

A xp(A a,I*j){
I d[a->r]; I w[a->r];
DO(a->r, d[i]=a->d[j[i]];
w[i]=a->w[j[i]])
A z=caa(a->p,r,d);
memcpy(z->w,w,a->r*sizeof(I)); R z; } /* transpose */

A sl0(A a,I x){ I r=a->r-1; A z=ah(r,a->d+1,0);
memcpy(z->w,a->w+1,r*sizeof(I));
z->p=a->p+x*a->w[0]; R z; }

A sl(A a,I*s,I*f){ I r=0;
DO(a->r, r+=s[i]!=f[i])
I d[r]; I w[r]; I j=0;
DO(r, while(s[j]==f[j])++j;
d[i]=1+(s[j]<f[j]? f[j]-s[j] : s[j]-f[j]);
w[i]= s[j]<f[j]? a->w[j] : -a->w[j];
++j)
A z=caa(a->p,r,d);
memcpy(z->w,w,r*sizeof(I)); DO(a->r, z->p += s[i] * a->w[i]) R z; } /* slices [s[i]..f[i]] */

I contig(A a){ I x=1;
DO(a->r, if(a->w[a->r-1-i]!=x) R 0;
x*=a->d[a->r-1-i])
R 1; }

pr(A a,I wid){ I max; I copy=0;
if (wid){ max=wid; } else {
I sz=tr(a->r,a->d);
if (!contig(a)){
a=cp(a);
copy=1;
}
max=0;
DO(sz,
sz=snprintf(NULL,0,"%d",a->p[i]);
if (sz>max) max=sz)
}
switch(a->r){
case 0: printf("%*d\n", max, *a->p); break;
case 1: {
I sep=0;
DO(a->d[0],
if (sep) printf(" ");
printf("%*d", max, *el(a,i));
sep=1) }
if (wid==0) printf("\n"); break;
default:DO(a->d[0], A t=sl0(a,i);
pr(t,max);
printf("\n");
free(t)) break;
}
if(copy)free(a); }

dmp(A a){
printf("%d\n",a->r);
DO(a->r,printf("%d ",a->d[i])) printf("\n");
DO(a->r,printf("%d ",a->w[i])) printf("\n");
printf("%p %p %p %p\n",
(void*)a, (void*)a->d, (void*)a->w, (void*)a->p); }

int main(){
A x=ar(3);
*el(x,0)=5;
*el(x,1)=6;
*el(x,2)=7;
dmp(x);
pr(x,0);

A a=iot(24);
dmp(a);
pr(a,0);
A b=ca(a->p,2,3,4);
dmp(b);
pr(b,0);
A c=sl(b,(I[]){0,1,1},(I[]){0,2,3});
pr(c,0);
A d=xp(c,(I[]){1,0});
pr(d,0);
return 0;
}

[luser droog]

On Monday, June 1, 2015 at 8:04:25 PM UTC-5, luser droog wrote:
> On Monday, June 1, 2015 at 3:36:14 AM UTC-5, Ian Collins wrote:
> > asetof...@gmail.com wrote:
> > > I wrote:"
> >
> > More gobbledegook?
> >
>
> If you like gobbledegook, here's my code all golfed-
> up and gook.

And buggy.
<snip>

> A xp(A a,I*j){
> I d[a->r]; I w[a->r];
> DO(a->r, d[i]=a->d[j[i]];
> w[i]=a->w[j[i]])
> A z=caa(a->p,r,d);

s.b. A z=caa(a->p,a->r,d);

[Chris M. Thomasson]

> > "luser droog" wrote in message
> > news:7beb4e90-e514-45e4...@googlegroups.com...

> On Monday, June 1, 2015 at 3:36:14 AM UTC-5, Ian Collins wrote:
> More gobbledegook?

> > If you like gobbledegook, here's my code all golfed-
> > up and gook. There's even a little whitespace in there.
> > A very little.

[...]

Holy Shi%!

;^o

[Rosario19]

What DO the macro above?
int i;
for(i=0; i<r; ++i)
{

while(s[j]==f[j])++j;
d[i]=1+(s[j]<f[j]? f[j]-s[j] : s[j]-f[j]);
w[i]= s[j]<f[j]? a->w[j] : -a->w[j];
++j

^^^^^^^^^^^^^^^^^^^^^
here at last there is need ';'
no it is ok... in macro DO() the ';' is in the macro
}

what is sl?

A new language?
I agree if operations ca, xp, sl, dmp, pr have a full names

[luser droog]

On Tuesday, June 2, 2015 at 1:30:21 AM UTC-5, Rosario19 wrote:
> On Mon, 1 Jun 2015 18:04:07 -0700 (PDT), luser droog
> <luser...@gmail.com> wrote:
>
> >On Monday, June 1, 2015 at 3:36:14 AM UTC-5, Ian Collins wrote:
> >> asetof...@gmail.com wrote:
> >> > I wrote:"
> >>
> >> More gobbledegook?
> >>
> >
> >If you like gobbledegook, here's my code all golfed-
> >up and gook. There's even a little whitespace in there.
> >A very little.
> >
> >$ cat arrg.c

[snip]

> >#define DO(n,x){I i=0,N=(n);for(;i<N;++i){x;}}/*loop i=[0..n)*/

[snip]

> >
> >A sl(A a,I*s,I*f){ I r=0;
> > DO(a->r, r+=s[i]!=f[i])
> > I d[r]; I w[r]; I j=0;
> > DO(r, while(s[j]==f[j])++j;
> > d[i]=1+(s[j]<f[j]? f[j]-s[j] : s[j]-f[j]);
> > w[i]= s[j]<f[j]? a->w[j] : -a->w[j];
> > ++j)
>
> What DO the macro above?
> int i;
> for(i=0; i<r; ++i)
> {
> while(s[j]==f[j])++j;
> d[i]=1+(s[j]<f[j]? f[j]-s[j] : s[j]-f[j]);
> w[i]= s[j]<f[j]? a->w[j] : -a->w[j];
> ++j
> ^^^^^^^^^^^^^^^^^^^^^
> here at last there is need ';'
> no it is ok... in macro DO() the ';' is in the macro
> }
>
> what is sl?

It takes slices from the array, returning a subarray.

[snip]

> >int main(){
> > A x=ar(3);
> > *el(x,0)=5;
> > *el(x,1)=6;
> > *el(x,2)=7;
> > dmp(x);
> > pr(x,0);
> >
> > A a=iot(24);
> > dmp(a);
> > pr(a,0);
> > A b=ca(a->p,2,3,4);
> > dmp(b);
> > pr(b,0);
> > A c=sl(b,(I[]){0,1,1},(I[]){0,2,3});
> > pr(c,0);
> > A d=xp(c,(I[]){1,0});
> > pr(d,0);
> > return 0;
> >}
> >
>
> A new language?
> I agree if operations ca, xp, sl, dmp, pr have a full names

It's a prototype for using dope vectors to manipulate array data.
It's related to my programming language 'inca'.
https://github.com/luser-dr00g/inca
But I'm not sure whether to incorporate this code into inca3.c
or start over and build on top of this.

There is a "readable" version here:
https://github.com/luser-dr00g/inca/blob/master/arrind.c
And an explanation on SO:
http://stackoverflow.com/questions/30409991/use-a-dope-vector-to-access-arbitrary-axial-slices-of-a-multidimensional-array