mirror of git://git.psyced.org/git/psyced
288 lines
6.8 KiB
C
288 lines
6.8 KiB
C
// $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;
|
|
}
|
|
|