eigenmath/simplify.cpp

387 lines
5.1 KiB
C++
Raw Permalink Normal View History

2004-05-28 20:08:11 +02:00
#include "stdafx.h"
2004-03-03 21:24:06 +01:00
#include "defs.h"
2008-05-18 21:34:03 +02:00
2004-05-28 20:08:11 +02:00
extern int trigmode;
2004-03-03 21:24:06 +01:00
static void simplify_tensor(void);
static int count(U *);
2005-06-25 21:29:07 +02:00
static int nterms(U *);
2004-05-28 20:08:11 +02:00
static void f1(void);
static void f2(void);
static void f3(void);
static void f4(void);
static void f5(void);
static void f9(void);
2004-03-03 21:24:06 +01:00
void
eval_simplify(void)
{
push(cadr(p1));
eval();
simplify();
}
void
simplify(void)
{
save();
2008-05-18 21:34:03 +02:00
simplify_main();
2005-07-22 01:23:51 +02:00
restore();
}
2008-05-18 21:34:03 +02:00
void
simplify_main(void)
2005-07-22 01:23:51 +02:00
{
2004-03-03 21:24:06 +01:00
p1 = pop();
2005-07-22 01:23:51 +02:00
if (istensor(p1)) {
2004-03-03 21:24:06 +01:00
simplify_tensor();
2005-07-22 01:23:51 +02:00
return;
}
if (find(p1, symbol(FACTORIAL))) {
2004-03-03 21:24:06 +01:00
push(p1);
2005-07-22 01:23:51 +02:00
simfac();
2006-04-19 02:51:30 +02:00
p2 = pop();
push(p1);
rationalize();
simfac();
p3 = pop();
if (count(p2) < count(p3))
p1 = p2;
else
p1 = p3;
2004-05-28 20:08:11 +02:00
}
2005-07-22 01:23:51 +02:00
f1();
f2();
f3();
f4();
f5();
f9();
push(p1);
2004-03-03 21:24:06 +01:00
}
static void
simplify_tensor(void)
{
int i;
p2 = alloc_tensor(p1->u.tensor->nelem);
p2->u.tensor->ndim = p1->u.tensor->ndim;
for (i = 0; i < p1->u.tensor->ndim; i++)
p2->u.tensor->dim[i] = p1->u.tensor->dim[i];
for (i = 0; i < p1->u.tensor->nelem; i++) {
push(p1->u.tensor->elem[i]);
simplify();
p2->u.tensor->elem[i] = pop();
}
if (iszero(p2))
p2 = zero; // null tensor becomes scalar zero
2004-03-03 21:24:06 +01:00
push(p2);
}
static int
count(U *p)
{
int n;
if (iscons(p)) {
n = 0;
while (iscons(p)) {
n += count(car(p)) + 1;
p = cdr(p);
}
} else
n = 1;
return n;
}
2004-05-28 20:08:11 +02:00
// try rationalizing
2004-03-03 21:24:06 +01:00
static void
2004-05-28 20:08:11 +02:00
f1(void)
2004-03-03 21:24:06 +01:00
{
2005-08-06 22:57:37 +02:00
if (car(p1) != symbol(ADD))
2004-03-03 21:24:06 +01:00
return;
push(p1);
rationalize();
p2 = pop();
if (count(p2) < count(p1))
p1 = p2;
}
2004-05-28 20:08:11 +02:00
// try condensing
2004-03-03 21:24:06 +01:00
static void
2004-05-28 20:08:11 +02:00
f2(void)
2004-03-03 21:24:06 +01:00
{
2005-08-06 22:57:37 +02:00
if (car(p1) != symbol(ADD))
2004-03-03 21:24:06 +01:00
return;
push(p1);
2006-01-15 00:55:11 +01:00
Condense();
2004-03-03 21:24:06 +01:00
p2 = pop();
if (count(p2) <= count(p1))
p1 = p2;
}
2004-05-28 20:08:11 +02:00
// this simplifies forms like (A-B) / (B-A)
2004-03-03 21:24:06 +01:00
static void
2004-05-28 20:08:11 +02:00
f3(void)
2004-03-03 21:24:06 +01:00
{
2004-05-28 20:08:11 +02:00
push(p1);
rationalize();
negate();
rationalize();
negate();
rationalize();
p2 = pop();
if (count(p2) < count(p1))
p1 = p2;
}
2004-03-03 21:24:06 +01:00
2004-05-28 20:08:11 +02:00
// try expanding denominators
2004-03-03 21:24:06 +01:00
2004-05-28 20:08:11 +02:00
static void
f4(void)
{
if (iszero(p1))
2004-03-03 21:24:06 +01:00
return;
2004-05-28 20:08:11 +02:00
push(p1);
rationalize();
inverse();
rationalize();
inverse();
rationalize();
2004-03-03 21:24:06 +01:00
p2 = pop();
if (count(p2) < count(p1))
p1 = p2;
}
2004-05-28 20:08:11 +02:00
// simplifies trig forms
2004-03-03 21:24:06 +01:00
2005-10-08 04:20:34 +02:00
void
simplify_trig(void)
{
save();
p1 = pop();
f5();
push(p1);
restore();
}
2004-03-03 21:24:06 +01:00
static void
2004-05-28 20:08:11 +02:00
f5(void)
2004-03-03 21:24:06 +01:00
{
2005-06-25 21:29:07 +02:00
if (find(p1, symbol(SIN)) == 0 && find(p1, symbol(COS)) == 0)
2004-03-03 21:24:06 +01:00
return;
p2 = p1;
trigmode = 1;
push(p2);
eval();
p3 = pop();
trigmode = 2;
push(p2);
eval();
p4 = pop();
trigmode = 0;
2005-06-25 21:29:07 +02:00
if (count(p4) < count(p3) || nterms(p4) < nterms(p3))
2004-03-03 21:24:06 +01:00
p3 = p4;
2005-06-25 21:29:07 +02:00
if (count(p3) < count(p1) || nterms(p3) < nterms(p1))
2004-03-03 21:24:06 +01:00
p1 = p3;
}
2004-05-28 20:08:11 +02:00
// if it's a sum then try to simplify each term
static void
f9(void)
{
if (car(p1) != symbol(ADD))
return;
push_integer(0);
p2 = cdr(p1);
while (iscons(p2)) {
push(car(p2));
simplify();
add();
p2 = cdr(p2);
}
p2 = pop();
if (count(p2) < count(p1))
p1 = p2;
}
2005-06-25 21:29:07 +02:00
static int
nterms(U *p)
{
if (car(p) != symbol(ADD))
return 1;
else
return length(p) - 1;
}
2007-05-08 16:57:30 +02:00
#if SELFTEST
2004-03-03 21:24:06 +01:00
static char *s[] = {
"simplify(A)",
"A",
"simplify(A+B)",
"A+B",
"simplify(A B)",
"A*B",
"simplify(A^B)",
"A^B",
"simplify(A/(A+B)+B/(A+B))",
"1",
"simplify((A-B)/(B-A))",
"-1",
"A=((A11,A12),(A21,A22))",
"",
"simplify(det(A) inv(A) - adj(A))",
"0",
"A=quote(A)",
"",
// this shows need for <= in try_factoring
2006-02-10 17:21:52 +01:00
// "x*(1+a)",
// "x+a*x",
2004-03-03 21:24:06 +01:00
2006-02-10 17:21:52 +01:00
// "simplify(last)",
// "x*(1+a)",
2004-03-03 21:24:06 +01:00
"simplify(-3 exp(-1/3 r + i phi) cos(theta) / sin(theta)"
" + 3 exp(-1/3 r + i phi) cos(theta) sin(theta)"
" + 3 exp(-1/3 r + i phi) cos(theta)^3 / sin(theta))",
"0",
2004-05-28 20:08:11 +02:00
"simplify((A^2 C^2 + A^2 D^2 + B^2 C^2 + B^2 D^2)/(A^2+B^2)/(C^2+D^2))",
"1",
2004-06-22 03:22:52 +02:00
"simplify(d(arctan(y/x),y))",
"x/(x^2+y^2)",
"simplify(d(arctan(y/x),x))",
"-y/(x^2+y^2)",
2005-06-25 21:29:07 +02:00
"simplify(1-sin(x)^2)",
"cos(x)^2",
"simplify(1-cos(x)^2)",
"sin(x)^2",
"simplify(sin(x)^2-1)",
"-cos(x)^2",
"simplify(cos(x)^2-1)",
"-sin(x)^2",
2006-04-19 02:51:30 +02:00
/*
2005-07-25 20:13:22 +02:00
"simfac(n!/n)-(n-1)!",
2005-07-22 00:41:09 +02:00
"0",
2005-07-25 20:13:22 +02:00
"simfac(n/n!)-1/(n-1)!",
2005-07-22 00:41:09 +02:00
"0",
2005-07-22 01:23:51 +02:00
2005-07-25 20:13:22 +02:00
"simfac(rationalize((n+k+1)/(n+k+1)!))-1/(n+k)!",
2005-07-22 01:23:51 +02:00
"0",
2005-07-25 20:13:22 +02:00
"simfac(condense((n+1)*n!))-(n+1)!",
2005-07-22 02:11:02 +02:00
"0",
2005-07-25 20:13:22 +02:00
"simfac(1/((n+1)*n!))-1/(n+1)!",
2005-07-22 01:23:51 +02:00
"0",
2005-07-22 02:22:44 +02:00
2005-07-25 20:13:22 +02:00
"simfac((n+1)!/n!)-n-1",
2005-07-22 02:22:44 +02:00
"0",
2005-07-25 20:13:22 +02:00
"simfac(n!/(n+1)!)-1/(n+1)",
2005-07-22 02:22:44 +02:00
"0",
2005-07-22 02:39:11 +02:00
2005-07-25 20:13:22 +02:00
"simfac(binomial(n+1,k)/binomial(n,k))",
2005-07-22 02:39:11 +02:00
"(1+n)/(1-k+n)",
2005-07-25 20:13:22 +02:00
"simfac(binomial(n,k)/binomial(n+1,k))",
2005-07-22 02:39:11 +02:00
"(1-k+n)/(1+n)",
2005-07-22 23:48:02 +02:00
2008-06-13 06:26:19 +02:00
"F(nn,kk)=kk*binomial(nn,kk)",
2005-07-22 23:48:02 +02:00
"",
2005-07-25 20:13:22 +02:00
"simplify(simfac((F(n,k)+F(n,k-1))/F(n+1,k))-n/(n+1))",
2005-07-22 23:48:02 +02:00
"0",
"F=quote(F)",
"",
2006-04-19 02:51:30 +02:00
*/
"simplify(n!/n)-(n-1)!",
"0",
"simplify(n/n!)-1/(n-1)!",
"0",
"simplify(rationalize((n+k+1)/(n+k+1)!))-1/(n+k)!",
"0",
"simplify(condense((n+1)*n!))-(n+1)!",
"0",
"simplify(1/((n+1)*n!))-1/(n+1)!",
"0",
"simplify((n+1)!/n!)-n-1",
"0",
"simplify(n!/(n+1)!)-1/(n+1)",
"0",
"simplify(binomial(n+1,k)/binomial(n,k))",
"(1+n)/(1-k+n)",
"simplify(binomial(n,k)/binomial(n+1,k))",
"(1-k+n)/(1+n)",
2008-06-13 06:26:19 +02:00
"F(nn,kk)=kk*binomial(nn,kk)",
2006-04-19 02:51:30 +02:00
"",
"simplify((F(n,k)+F(n,k-1))/F(n+1,k))-n/(n+1)",
"0",
"F=quote(F)",
"",
"simplify((n+1)/(n+1)!)-1/n!",
"0",
2006-09-19 23:04:36 +02:00
"simplify(a*b+a*c)",
"a*(b+c)",
2008-05-18 21:34:03 +02:00
// Symbol's binding is evaluated, undoing simplify
2006-09-19 23:04:36 +02:00
"x=simplify(a*b+a*c)",
"",
"x",
2008-05-18 21:34:03 +02:00
"a*b+a*c",
2004-03-03 21:24:06 +01:00
};
void
test_simplify(void)
{
test(__FILE__, s, sizeof s / sizeof (char *));
}
2007-05-08 16:57:30 +02:00
#endif