123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298 |
- #include "stdafx.h"
- #include "defs.h"
- static int
- check_arg(void)
- {
- if (!istensor(p1))
- return 0;
- else if (p1->u.tensor->ndim != 2)
- return 0;
- else if (p1->u.tensor->dim[0] != p1->u.tensor->dim[1])
- return 0;
- else
- return 1;
- }
- void
- det(void)
- {
- int i, n;
- U **a;
- save();
- p1 = pop();
- if (check_arg() == 0) {
- push_symbol(DET);
- push(p1);
- list(2);
- restore();
- return;
- }
- n = p1->u.tensor->nelem;
- a = p1->u.tensor->elem;
- for (i = 0; i < n; i++)
- if (!isnum(a[i]))
- break;
- if (i == n)
- yydetg();
- else {
- for (i = 0; i < p1->u.tensor->nelem; i++)
- push(p1->u.tensor->elem[i]);
- determinant(p1->u.tensor->dim[0]);
- }
- restore();
- }
- void
- determinant(int n)
- {
- int h, i, j, k, q, s, sign, t;
- int *a, *c, *d;
- h = tos - n * n;
- a = (int *) malloc(3 * n * sizeof (int));
- if (a == NULL)
- out_of_memory();
- c = a + n;
- d = c + n;
- for (i = 0; i < n; i++) {
- a[i] = i;
- c[i] = 0;
- d[i] = 1;
- }
- sign = 1;
- push(zero);
- for (;;) {
- if (sign == 1)
- push_integer(1);
- else
- push_integer(-1);
- for (i = 0; i < n; i++) {
- k = n * a[i] + i;
- push(stack[h + k]);
- multiply();
- }
- add();
-
- j = n - 1;
- s = 0;
- P4: q = c[j] + d[j];
- if (q < 0) {
- d[j] = -d[j];
- j--;
- goto P4;
- }
- if (q == j + 1) {
- if (j == 0)
- break;
- s++;
- d[j] = -d[j];
- j--;
- goto P4;
- }
- t = a[j - c[j] + s];
- a[j - c[j] + s] = a[j - q + s];
- a[j - q + s] = t;
- c[j] = q;
- sign = -sign;
- }
- free(a);
- stack[h] = stack[tos - 1];
- tos = h + 1;
- }
- void
- detg(void)
- {
- save();
- p1 = pop();
- if (check_arg() == 0) {
- push_symbol(DET);
- push(p1);
- list(2);
- restore();
- return;
- }
- yydetg();
- restore();
- }
- void
- yydetg(void)
- {
- int i, n;
- n = p1->u.tensor->dim[0];
- for (i = 0; i < n * n; i++)
- push(p1->u.tensor->elem[i]);
- lu_decomp(n);
- tos -= n * n;
- push(p1);
- }
- #define M(i, j) stack[h + n * (i) + (j)]
- void
- lu_decomp(int n)
- {
- int d, h, i, j;
- h = tos - n * n;
- p1 = one;
- for (d = 0; d < n - 1; d++) {
-
- if (equal(M(d, d), zero)) {
-
- for (i = d + 1; i < n; i++)
- if (!equal(M(i, d), zero))
- break;
- if (i == n) {
- p1 = zero;
- break;
- }
-
- for (j = d; j < n; j++) {
- p2 = M(d, j);
- M(d, j) = M(i, j);
- M(i, j) = p2;
- }
-
- push(p1);
- negate();
- p1 = pop();
- }
-
- push(p1);
- push(M(d, d));
- multiply();
- p1 = pop();
-
- for (i = d + 1; i < n; i++) {
-
- push(M(i, d));
- push(M(d, d));
- divide();
- negate();
- p2 = pop();
-
- M(i, d) = zero;
- for (j = d + 1; j < n; j++) {
- push(M(d, j));
- push(p2);
- multiply();
- push(M(i, j));
- add();
- M(i, j) = pop();
- }
- }
- }
-
- push(p1);
- push(M(n - 1, n - 1));
- multiply();
- p1 = pop();
- }
|