202 lines
2.4 KiB
C++
202 lines
2.4 KiB
C++
#include "stdafx.h"
|
|
|
|
|
|
|
|
#include "defs.h"
|
|
|
|
extern int powermode;
|
|
|
|
static int check_args(void);
|
|
static void charpoly_n(void);
|
|
static void charpoly_s(void);
|
|
|
|
#define A p1
|
|
#define B p2
|
|
#define LAMBDA p3
|
|
#define C p4
|
|
#define POLY p5
|
|
|
|
void
|
|
charpoly(void)
|
|
{
|
|
int i, n;
|
|
U **a;
|
|
save();
|
|
LAMBDA = pop();
|
|
A = pop();
|
|
if (check_args() == 0) {
|
|
push_symbol(CHARPOLY);
|
|
push(A);
|
|
push(LAMBDA);
|
|
list(3);
|
|
restore();
|
|
return;
|
|
} else {
|
|
n = A->u.tensor->nelem;
|
|
a = A->u.tensor->elem;
|
|
for (i = 0; i < n; i++)
|
|
if (!isnum(a[i]))
|
|
break;
|
|
if (i == n)
|
|
charpoly_n();
|
|
else
|
|
charpoly_s();
|
|
}
|
|
|
|
// multiply constant terms by LAMBDA ^ 0
|
|
|
|
push_symbol(POWER);
|
|
push(LAMBDA);
|
|
push(_zero);
|
|
list(3);
|
|
multiply();
|
|
|
|
restore();
|
|
}
|
|
|
|
static int
|
|
check_args(void)
|
|
{
|
|
if (!istensor(A))
|
|
return 0;
|
|
else if (A->u.tensor->ndim != 2)
|
|
return 0;
|
|
else if (A->u.tensor->dim[0] != A->u.tensor->dim[1])
|
|
return 0;
|
|
else if (!issym(LAMBDA))
|
|
return 0;
|
|
else
|
|
return 1;
|
|
}
|
|
|
|
static void
|
|
charpoly_n(void)
|
|
{
|
|
int i, j, n;
|
|
|
|
n = A->u.tensor->dim[0];
|
|
|
|
// copy A
|
|
|
|
push(A);
|
|
eval();
|
|
B = pop();
|
|
|
|
// C = trace(A)
|
|
|
|
push(A);
|
|
trace();
|
|
C = pop();
|
|
|
|
// POLY = LAMBDA ^ n - C * LAMBDA ^ (n - 1)
|
|
|
|
push(LAMBDA);
|
|
push_integer(n);
|
|
power();
|
|
push(C);
|
|
push(LAMBDA);
|
|
push_integer(n - 1);
|
|
power();
|
|
multiply();
|
|
subtract();
|
|
POLY = pop();
|
|
|
|
for (i = 1; i < n; i++) {
|
|
|
|
// B = B - C * I
|
|
|
|
for (j = 0; j < n; j++) {
|
|
push(B->u.tensor->elem[j * n + j]);
|
|
push(C);
|
|
subtract();
|
|
B->u.tensor->elem[j * n + j] = pop();
|
|
}
|
|
|
|
// B = A . B
|
|
|
|
push(A);
|
|
push(B);
|
|
inner();
|
|
B = pop();
|
|
|
|
// C = trace(B) / (i + 1)
|
|
|
|
push(B);
|
|
trace();
|
|
push_integer(i + 1);
|
|
divide();
|
|
C = pop();
|
|
|
|
// POLY = POLY - C * LAMBDA ^ (n - i - 1)
|
|
|
|
push(POLY);
|
|
push(C);
|
|
push(LAMBDA);
|
|
push_integer(n - i - 1);
|
|
power();
|
|
multiply();
|
|
subtract();
|
|
POLY = pop();
|
|
}
|
|
|
|
push(POLY);
|
|
|
|
if (n % 2)
|
|
negate();
|
|
}
|
|
|
|
static void
|
|
charpoly_s(void)
|
|
{
|
|
int n;
|
|
n = A->u.tensor->dim[0];
|
|
push(A);
|
|
push(LAMBDA);
|
|
push_identity_matrix(n);
|
|
multiply();
|
|
subtract();
|
|
det();
|
|
}
|
|
|
|
static char *s[] = {
|
|
|
|
"A=((1,2),(3,4))",
|
|
"",
|
|
|
|
"f(x)=eval(charpoly(A))",
|
|
"",
|
|
|
|
"f",
|
|
"-5*x-2*x^0+x^2",
|
|
|
|
"f(A)",
|
|
"0",
|
|
|
|
"A=((a,b),(c,d))",
|
|
"",
|
|
|
|
"f(lambda)=eval(charpoly(A,lambda))",
|
|
"",
|
|
|
|
"f",
|
|
"a*d*lambda^0-a*lambda-b*c*lambda^0-d*lambda+lambda^2",
|
|
|
|
"f(A)",
|
|
"0",
|
|
|
|
"A=quote(A)",
|
|
"",
|
|
|
|
"f()=eval(quote(f))",
|
|
"",
|
|
|
|
"charpoly(A)",
|
|
"charpoly(A,x)",
|
|
};
|
|
|
|
void
|
|
test_charpoly(void)
|
|
{
|
|
test(__FILE__, s, sizeof s / sizeof (char *));
|
|
}
|