#pragma	src	"/usr/doc/src/linalg"
#pragma	lib	"liblinalg.a"

typedef float real;
extern const float meps;

/* Cubes, matrics, and Vector data types */

/* Morton ordering */
#define evenbits 0xAAAAAAAA
#define oddbits 0x55555555
//#define threebits 0x49249249
#define oinc2(i) (((i) - oddbits) & oddbits)
#define einc2(i) (((i) - evenbits) & evenbits)
#define odec2(i) (((i) - 1) & oddbits)
#define edec2(i) (((i) - 1) & evenbits)
#define tr(i, j) (((i) << 1) | ((j) >> 1))
ushort undilate2(uint);
uint dilate2(ushort);

/* Iterators */
#define miter(m, i) for((i) = 0; (i) < (m)->len; (i)++) if(mbndok((m), (i)))
#define viter(v, i) for((i) = 0; (i) < (v)->len; (i)++)

typedef struct block block;
typedef struct cube cube;
typedef struct matrix matrix;
typedef struct vector vector;

#define C(x, i, j, k) ((x)->data[(k) * (x)->stride + (dilate2(i) | (dilate2(j) << 1))])
#define CV(x, i, k) ((x)->data[(k) * (x)->stride + (i)])
#define M(x, i, j) ((x)->data[dilate2(i) | (dilate2(j) << 1)])
#define V(x, i) ((x)->data[(i)])

#define mbndok(a, i) (((i) & oddbits) < (a)->mr && ((i) & evenbits) < (a)->mc)

struct block{
	real *data;
	ulong sz;
	uchar allocated:1;
};

struct vector {
	block;
	uint len;
};

struct matrix{
	vector;
	ushort nr, nc;
	uint mr, mc;
	uint levels, offset;
};

struct cube{
	matrix;
	ushort ns;
	uint stride;
};

cube *ncube(cube *, real *, ushort, ushort, ushort);
matrix *nmatrix(matrix *, real *, ushort, ushort);
vector *nvector(vector *, real *, uint);

void fblock(block *);
block *zblock(block *);
block *cblock(block *, block *);

matrix *qdot(matrix *, matrix *, matrix *, int);
matrix *mdot(matrix *, matrix *, matrix *, int);
matrix *madd(matrix *, matrix *, matrix *);
matrix *smprod(real, matrix *, matrix *);
real dot(vector *, vector *);
vector *sprod(real, vector *, vector *);
vector *vadd(vector *, vector *, vector *);
vector *vprod(vector *, vector *, vector *);
vector *vsub(vector *, vector *, vector *);
vector *vdot(matrix *, vector *, vector *, int);
void vnorm(uint, ...);

void vsort(vector *);
uint *ivsort(vector *);

extern uint MAXITERS;
matrix *eigen(matrix *, matrix *, vector *);
matrix *jacobi(matrix *, matrix *, vector *);
matrix *tridiag(matrix *, matrix *, matrix *);
matrix *nystrom(matrix *, matrix *, vector *);
matrix *qr(matrix *, matrix *, matrix *);
matrix *gs(matrix *);

real entropy(matrix *);
matrix *mahalanobis(matrix *, matrix *);
matrix *msqrt(matrix *, matrix *);
matrix *minv(matrix *, matrix *, real);
vector *solve(matrix *, vector *, vector *);

void svd(matrix *, matrix *, vector *, matrix *);

matrix *col(matrix *, matrix *, uint);
matrix *row(matrix *, matrix *, uint);
matrix *slice(cube *, matrix *, uint);
vector *drop(matrix *, vector *);
matrix *raise(vector *, matrix *);
matrix *trans(matrix *, matrix *);
matrix *diag(vector *, matrix *);
vector *msum(matrix *, vector *, uint);
vector *getdiag(matrix *, vector *);