psyced/world/net/matrix.c

288 lines
6.8 KiB
C
Raw Permalink Normal View History

// $Id: matrix.c,v 1.10 2005/03/14 10:23:26 lynx Exp $ // vim:syntax=lpc
//
// just a joke, really
//
#include <net.h>
#define CL_INDEX (#'[)
#define CL_IF (#'?)
#define CL_NIF (#'?!)
#define CL_RANGE (#'[..])
#define CL_L_RANGE (#'[..)
protected int *dim(mixed M);
protected mixed TRANS(mixed M);
protected varargs mixed MpM(mixed A, mixed B, mixed C);
protected varargs mixed MxM(mixed A, mixed B, mixed C);
// delta(ij)
protected mixed E(int n);
// multipliziere die i-te Zeile mit x.. Ri(x) element M(nxn,K)
protected mixed R(int n, int i, int x);
// addiere das xfache der i-ten Zeile zur j-ten Zeile.. Qij(x) element M(nxn,K)
// WATCH OUT: von links!
protected mixed Q(int n, int i, int x, int j);
// vertauscht i-te und j-te zeile
protected mixed P(int n, int i, int j);
// BEWARE!! the mysterious sparcematrix optimizer! (M -> L)
//
// transforms a matrix into its aequivalent transformation.
// afterwards it can be used one-way.. funcall(L, some-matrix);
// the transformed matrix is much faster in case it contains many zeros
// and zero-lines.
//
//
//
protected closure matrix2closure(mixed matrix);
// reverse
protected mixed closure2matrix(closure c);
protected string printMatrix(mixed matrix);
protected int spur(mixed matrix);
/*
* Matrizen in der Form M[i][j].. , das ist per
* Definition Mij. Also
* i Zeilen
* j Spalten
*
*
*/
protected int *dim(mixed M) {
int *dim, *row;
if (closurep(M)) return funcall(M,"dim");
dim = ({ 0,0 });
dim[0] = sizeof(M);
dim[1] = sizeof(M[0]);
if(dim[0] == 0 || dim[1] == 0) return ({ -1,-1 });
foreach (row : M) {
if (sizeof(row) != dim[1])
return ({ -1,-1 });
}
return dim;
}
protected mixed TRANS(mixed M) {
mixed *T;
int *dimM, wantClosure, i, j;
if(closurep(M)) {
dimM = funcall(M,"dim");
M = closure2matrix(M);
wantClosure = 1;
} else dimM = dim(M);
if(dimM[0] == -1) return ({ -1,-1 });
T = allocate( ({ dimM[1],dimM[0] }) );
for (i = 0; i <= dimM[0] - 1; i++) { // loop rows of M
for (j = 0; j <= dimM[1] - 1; j++) { // loop cols of M
T[j][i] = M[i][j];
}
}
unless (wantClosure) return T;
return matrix2closure(T);
}
// can be used as A * B = C
protected varargs mixed MxM(mixed A, mixed B, mixed C) {
int *dimA, *dimB, *dimC, sum, i, j, k;
if(closurep(B)) {
dimB = funcall(B,"dim");
B = closure2matrix(B);
} else dimB = dim(B);
if(closurep(A)) {
return funcall(A,B);
}
dimA = dim(A); // get dimensions to check whether
if (pointerp(C)) {
dimC = dim(C);
if (dimC[0] != dimB[0] || dimC[1] != dimA[1] || dimC[0] == -1) return C;
} else {
C = allocate( ({ dimA[0], dimB[1] }) );
}
if (dimA[1] != dimB[0] || dimB[0] == -1 || dimA[0] == -1) return C;
// matrices are valid, we can start multiplication
for (i = 0; i <= dimA[0] - 1; i++) { // loop rows of A
for (j = 0; j <= dimB[1] - 1; j++) { // loop cols of B
sum = 0;
for (k = 0; k <= dimA[1] - 1; k++) { // vektor-produkt
sum += A[i][k] * B[k][j];
}
C[i][j] = sum;
}
}
return C;
}
protected varargs mixed MpM(mixed A, mixed B, mixed C) {
int *dimA, *dimB, *dimC, j, i;
if (closurep(A)) {
dimA = funcall(A,"dim");
A = closure2matrix(A);
} else dimA = dim(A);
if (closurep(B)) {
dimB = funcall(B,"dim");
B = closure2matrix(B);
} else dimB = dim(B);
if (pointerp(C)) {
dimC = dim(C);
if (dimC[0] != dimB[0] || dimC[1] != dimB[1]) return C;
} else {
C = allocate( ({ dimB[0], dimB[1] }) );
}
if (dimA[0] != dimB[0] || dimA[1] != dimB[1]) return C;
// matrices are valid, lets start add
for (i = 0; i <= dimA[0] - 1; i++) { // loop
for (j = 0; j <= dimB[1] - 1; j++) { // loop
C[i][j] = B[i][j] + A[i][j];
}
}
return C;
}
protected mixed E(int n) {
int i;
mixed matrix;
matrix = allocate(({ n,n }));
for(i = 0; i < n; i++) {
matrix[i][i] = 1;
}
return matrix;
}
protected mixed R(int n, int i, int x) {
mixed matrix;
i--;
matrix = E(n);
if (i < 0 || i >= n)
return matrix;
matrix[i][i] = x;
return matrix;
}
protected mixed Q(int n, int i, int x, int j) {
mixed matrix;
i--; j--;
matrix = E(n);
if (i == j || i < 0 || j < 0 || j >= n || i >= n)
return matrix;
// zeile j, spalte i.. Qij(x)
matrix[j][i] = x;
return matrix;
}
// exchanges j'th and i'th rows of the matrix
protected mixed P(int n, int i, int j) {
mixed matrix;
i--; j--;
matrix = E(n);
if (i == j || i < 0 || j < 0 || j >= n || i >= n)
return matrix;
matrix[i][i] = 0;
matrix[j][i] = 1;
matrix[i][j] = 1;
matrix[j][j] = 0;
return matrix;
}
// closure gets the matrix.. we are working on a symbol!
// 'matrix
// in column 'j
//
protected closure matrix2closure(mixed matrix) {
int *dim = dim(matrix);
int flag, i, j;
mixed *c_array = ({ (#',) });
mixed *temp;
for (i = 0; i <= dim[0] - 1; i++) {
flag = 0;
temp = 0;
for (j = 0; j <= dim[1] - 1; j++) {
switch (matrix[i][j]) {
case 0:
break;
case 1:
flag = 1;
temp = ({ (#'+), ({ CL_INDEX, ({ CL_INDEX, 'matrix, j }), 'j }), temp });
break;
default:
flag = 1;
temp = ({ (#'+), ({ #'*, matrix[i][j], ({ CL_INDEX, ({ CL_INDEX, 'matrix, j }), 'j }) }), temp });
break;
}
}
if (flag == 1) {
c_array += ({ ({ #'=, ({ CL_INDEX, ({ CL_INDEX, 'matrix_out, i }), 'j }), temp }) });
}
}
// P2(("%O",c_array))
return lambda(({ 'matrix }),
({ (#',),
({ CL_IF, ({ #'==, 'matrix, "dim" }), ({ #'return, quote(dim) }) }),
({ #'=, 'dim, ({ symbol_function("dim", ME), 'matrix }) }),
({ #'=, 'alloc_dim, ({ #'({, dim[0], ({ CL_INDEX, 'dim, 1 }) }) }),
// ({ #'=, ({ CL_INDEX, 'alloc_dim, 1 }), ({ CL_INDEX, 'dim, 1 }) }),
({ #'=, 'matrix_out, ({ #'allocate, 'alloc_dim }) }),
({ CL_IF, ({ #'||, ({ #'==, ({ CL_INDEX, 'dim, 0 }), -1 }),
({ #'!=, ({ CL_INDEX, 'dim, 0 }), dim[1] }) }),
({ #'return, 'matrix_out })
}),
({ #'=, 'j, 0}),
({ #'while,
({ #'<= , 'j, ({ #'-, ({ CL_INDEX, 'dim, 1 }), 1 }) }),
({ #'return, 'matrix_out }),
c_array,
({ #'++, 'j})
})
})
);
}
// just multiply with E
protected mixed closure2matrix(closure c) {
return funcall(c, E( funcall(c,"dim")[1] ) );
}
protected string printMatrix(mixed matrix) {
string output;
int *row;
if (closurep(matrix)) matrix = closure2matrix(matrix);
output = "";
foreach (row : matrix) {
output += "|\t"+ implode(map(row,#'to_string), "\t") +"\t|\n";
}
return output;
}
protected int spur(mixed matrix) {
int *dim, c, n, spur;
spur = 0;
dim = dim(matrix);
n = min(dim[0], dim[1]) - 1;
for (c = 0; c <= n; c++) {
spur += matrix[c][c];
}
return spur;
}