eigenmath/expand.cpp

609 lines
9.0 KiB
C++
Raw Permalink Normal View History

2008-05-18 18:31:15 +02:00
// Partial fraction expansion
2008-06-15 02:05:24 +02:00
//
2008-06-15 18:28:46 +02:00
// Example
//
// expand(1/(x^3+x^2),x)
2008-06-15 02:05:24 +02:00
//
// 1 1 1
// ---- - --- + -------
// 2 x x + 1
// x
2008-05-18 18:31:15 +02:00
2008-04-12 05:28:17 +02:00
#include "stdafx.h"
#include "defs.h"
void
eval_expand(void)
{
2008-06-17 04:33:14 +02:00
// 1st arg
2008-06-15 02:05:24 +02:00
push(cadr(p1));
2008-04-12 05:28:17 +02:00
eval();
2008-06-17 04:33:14 +02:00
// 2nd arg
2008-06-15 02:05:24 +02:00
push(caddr(p1));
2008-04-12 05:28:17 +02:00
eval();
2008-06-17 04:33:14 +02:00
2008-04-12 05:28:17 +02:00
p2 = pop();
if (p2 == symbol(NIL))
guess();
else
push(p2);
2008-06-17 04:33:14 +02:00
2008-04-12 05:28:17 +02:00
expand();
}
2008-04-14 02:20:00 +02:00
#define A p2
#define B p3
2008-06-07 04:54:52 +02:00
#define C p4
#define F p5
#define P p6
#define Q p7
#define T p8
#define X p9
2008-04-12 05:28:17 +02:00
void
expand(void)
{
save();
X = pop();
F = pop();
2008-04-19 03:28:05 +02:00
if (istensor(F)) {
expand_tensor();
restore();
return;
}
2008-04-12 05:28:17 +02:00
// if sum of terms then sum over the expansion of each term
if (car(F) == symbol(ADD)) {
push_integer(0);
p1 = cdr(F);
while (iscons(p1)) {
push(car(p1));
push(X);
expand();
add();
p1 = cdr(p1);
}
restore();
return;
}
2008-04-13 23:36:03 +02:00
// B = numerator
2008-04-12 05:28:17 +02:00
push(F);
numerator();
2008-04-13 23:36:03 +02:00
B = pop();
2008-06-07 04:54:52 +02:00
// A = denominator
2008-04-12 05:28:17 +02:00
push(F);
denominator();
2008-04-13 23:36:03 +02:00
A = pop();
2008-04-12 05:28:17 +02:00
2008-06-07 04:54:52 +02:00
remove_negative_exponents();
2008-04-19 04:24:34 +02:00
2008-06-07 04:54:52 +02:00
// Q = quotient
2008-04-12 05:28:17 +02:00
push(B);
2008-04-13 23:36:03 +02:00
push(A);
2008-04-12 05:28:17 +02:00
push(X);
divpoly();
Q = pop();
2008-06-07 04:54:52 +02:00
// remainder B = B - A * Q
2008-04-12 05:28:17 +02:00
push(B);
2008-04-13 23:36:03 +02:00
push(A);
2008-04-12 05:28:17 +02:00
push(Q);
multiply();
subtract();
2008-04-13 23:36:03 +02:00
B = pop();
2008-04-12 05:28:17 +02:00
// if the remainder is zero then we're done
2008-04-13 23:36:03 +02:00
if (iszero(B)) {
2008-04-14 04:35:01 +02:00
push(Q);
2008-04-12 05:28:17 +02:00
restore();
return;
}
2008-06-07 18:07:10 +02:00
// A = factor(A)
2008-04-12 05:28:17 +02:00
2008-06-07 04:54:52 +02:00
push(A);
push(X);
factorpoly();
A = pop();
2008-04-14 04:35:01 +02:00
2008-06-07 18:07:10 +02:00
expand_get_C();
expand_get_B();
expand_get_A();
2008-06-07 03:15:47 +02:00
2008-06-07 04:54:52 +02:00
if (istensor(C)) {
push(C);
inv();
push(B);
inner();
push(A);
inner();
} else {
push(B);
push(C);
divide();
push(A);
multiply();
}
2008-06-07 03:15:47 +02:00
push(Q);
add();
2008-04-12 05:28:17 +02:00
restore();
}
2008-04-19 03:28:05 +02:00
void
expand_tensor(void)
{
int i;
push(F);
copy_tensor();
F = pop();
for (i = 0; i < F->u.tensor->nelem; i++) {
push(F->u.tensor->elem[i]);
push(X);
expand();
F->u.tensor->elem[i] = pop();
}
push(F);
}
2008-04-19 04:24:34 +02:00
void
2008-06-07 04:54:52 +02:00
remove_negative_exponents(void)
2008-04-19 04:24:34 +02:00
{
2008-06-07 04:54:52 +02:00
int h, i, j, k, n;
2008-04-19 04:24:34 +02:00
2008-06-07 04:54:52 +02:00
h = tos;
factors(A);
factors(B);
n = tos - h;
2008-04-19 04:24:34 +02:00
// find the smallest exponent
2008-06-07 04:54:52 +02:00
j = 0;
for (i = 0; i < n; i++) {
p1 = stack[h + i];
2008-04-19 04:24:34 +02:00
if (car(p1) != symbol(POWER))
continue;
if (cadr(p1) != X)
continue;
push(caddr(p1));
2008-06-07 04:54:52 +02:00
k = pop_integer();
if (k == (int) 0x80000000)
2008-04-19 04:27:57 +02:00
continue;
2008-06-07 04:54:52 +02:00
if (k < j)
j = k;
2008-04-19 04:24:34 +02:00
}
2008-06-07 04:54:52 +02:00
tos = h;
2008-04-19 04:24:34 +02:00
2008-06-07 04:54:52 +02:00
if (j == 0)
2008-04-19 04:24:34 +02:00
return;
2008-06-07 05:27:43 +02:00
// A = A / X^j
2008-04-19 04:24:34 +02:00
push(A);
push(X);
2008-06-07 05:27:43 +02:00
push_integer(-j);
2008-04-19 04:24:34 +02:00
power();
multiply();
A = pop();
2008-04-14 04:35:01 +02:00
2008-06-07 05:27:43 +02:00
// B = B / X^j
2008-04-14 04:35:01 +02:00
2008-06-07 04:54:52 +02:00
push(B);
2008-04-14 04:35:01 +02:00
push(X);
2008-06-07 05:27:43 +02:00
push_integer(-j);
2008-06-07 04:54:52 +02:00
power();
multiply();
B = pop();
2008-04-14 04:35:01 +02:00
}
2008-06-07 04:54:52 +02:00
// Returns the expansion coefficient matrix C.
2008-06-07 18:07:10 +02:00
//
2008-06-07 18:14:39 +02:00
// Example:
2008-06-07 18:07:10 +02:00
//
// B 1
// --- = -----------
// A 2
// x (x + 1)
//
2008-06-07 18:14:39 +02:00
// We have
2008-06-07 18:07:10 +02:00
//
// B Y1 Y2 Y3
// --- = ---- + ---- + -------
// A 2 x x + 1
// x
//
2008-06-07 18:14:39 +02:00
// Our task is to solve for the unknowns Y1, Y2, and Y3.
//
// Multiplying both sides by A yields
2008-06-07 18:07:10 +02:00
//
// AY1 AY2 AY3
// B = ----- + ----- + -------
// 2 x x + 1
// x
//
// Let
//
// A A A
// W1 = ---- W2 = --- W3 = -------
// 2 x x + 1
// x
//
// Then the coefficient matrix C is
//
// coeff(W1,x,0) coeff(W2,x,0) coeff(W3,x,0)
//
// C = coeff(W1,x,1) coeff(W2,x,1) coeff(W3,x,1)
//
// coeff(W1,x,2) coeff(W2,x,2) coeff(W3,x,2)
2008-06-07 23:56:54 +02:00
//
// It follows that
//
// coeff(B,x,0) Y1
//
// coeff(B,x,1) = C Y2
//
// coeff(B,x,2) = Y3
//
// Hence
//
// Y1 coeff(B,x,0)
// -1
// Y2 = C coeff(B,x,1)
//
// Y3 coeff(B,x,2)
2008-04-14 02:20:00 +02:00
2008-06-07 04:54:52 +02:00
void
2008-06-07 18:07:10 +02:00
expand_get_C(void)
2008-04-12 05:28:17 +02:00
{
2008-06-07 04:54:52 +02:00
int h, i, j, n;
U **a;
h = tos;
2008-06-07 03:15:47 +02:00
if (car(A) == symbol(MULTIPLY)) {
p1 = cdr(A);
while (iscons(p1)) {
F = car(p1);
2008-06-07 18:07:10 +02:00
expand_get_CF();
2008-06-07 03:15:47 +02:00
p1 = cdr(p1);
}
} else {
F = A;
2008-06-07 18:07:10 +02:00
expand_get_CF();
2008-06-07 03:15:47 +02:00
}
2008-06-07 04:54:52 +02:00
n = tos - h;
if (n == 1) {
C = pop();
return;
}
C = alloc_tensor(n * n);
C->u.tensor->ndim = 2;
C->u.tensor->dim[0] = n;
C->u.tensor->dim[1] = n;
a = stack + h;
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
push(a[j]);
push(X);
push_integer(i);
power();
divide();
push(X);
filter();
C->u.tensor->elem[n * i + j] = pop();
}
}
tos -= n;
2008-06-07 03:15:47 +02:00
}
2008-04-13 23:36:03 +02:00
2008-06-07 18:07:10 +02:00
// The following table shows the push order for simple roots, repeated roots,
// and inrreducible factors.
//
// Factor F Push 1st Push 2nd Push 3rd Push 4th
//
//
// A
// x ---
// x
2008-06-07 03:15:47 +02:00
//
//
2008-06-07 18:07:10 +02:00
// 2 A A
// x ---- ---
// 2 x
// x
2008-06-07 03:15:47 +02:00
//
//
2008-06-07 18:07:10 +02:00
// A
// x + 1 -------
// x + 1
2008-06-07 03:15:47 +02:00
//
//
2008-06-07 18:07:10 +02:00
// 2 A A
// (x + 1) ---------- -------
// 2 x + 1
// (x + 1)
2008-06-07 03:15:47 +02:00
//
//
2008-06-07 18:07:10 +02:00
// 2 A Ax
// x + x + 1 ------------ ------------
// 2 2
// x + x + 1 x + x + 1
2008-06-07 03:15:47 +02:00
//
//
2008-06-07 18:07:10 +02:00
// 2 2 A Ax A Ax
// (x + x + 1) --------------- --------------- ------------ ------------
// 2 2 2 2 2 2
// (x + x + 1) (x + x + 1) x + x + 1 x + x + 1
2008-06-07 03:15:47 +02:00
//
//
2008-06-07 18:07:10 +02:00
// For T = A/F and F = P^N we have
2008-06-07 03:15:47 +02:00
//
//
2008-06-07 18:07:10 +02:00
// Factor F Push 1st Push 2nd Push 3rd Push 4th
2008-06-07 03:15:47 +02:00
//
2008-06-07 18:07:10 +02:00
// x T
2008-06-07 03:15:47 +02:00
//
2008-06-07 18:07:10 +02:00
// 2
// x T TP
2008-06-07 03:15:47 +02:00
//
2008-06-07 23:56:54 +02:00
//
2008-06-07 18:07:10 +02:00
// x + 1 T
2008-06-07 03:15:47 +02:00
//
2008-06-07 18:07:10 +02:00
// 2
// (x + 1) T TP
2008-06-07 03:15:47 +02:00
//
2008-06-07 18:07:10 +02:00
// 2
// x + x + 1 T TX
//
// 2 2
// (x + x + 1) T TX TP TPX
//
//
// Hence we want to push in the order
//
// T * (P ^ i) * (X ^ j)
2008-06-07 03:15:47 +02:00
//
// for all i, j such that
//
2008-06-07 18:07:10 +02:00
// i = 0, 1, ..., N - 1
//
// j = 0, 1, ..., deg(P) - 1
2008-06-07 03:15:47 +02:00
//
2008-06-07 18:07:10 +02:00
// where index j runs first.
2008-04-13 23:36:03 +02:00
2008-06-07 03:15:47 +02:00
void
2008-06-07 18:07:10 +02:00
expand_get_CF(void)
2008-06-07 03:15:47 +02:00
{ int d, i, j, n;
if (!find(F, X))
return;
trivial_divide();
2008-06-07 18:07:10 +02:00
if (car(F) == symbol(POWER)) {
2008-06-07 03:15:47 +02:00
push(caddr(F));
n = pop_integer();
P = cadr(F);
} else {
n = 1;
P = F;
}
push(P);
push(X);
degree();
d = pop_integer();
for (i = 0; i < n; i++) {
for (j = 0; j < d; j++) {
push(T);
push(P);
push_integer(i);
power();
multiply();
push(X);
push_integer(j);
power();
multiply();
}
}
}
2008-04-13 23:36:03 +02:00
2008-06-07 05:27:43 +02:00
// Returns T = A/F where F is a factor of A.
2008-04-14 02:20:00 +02:00
2008-06-07 03:15:47 +02:00
void
trivial_divide(void)
{
int h;
if (car(A) == symbol(MULTIPLY)) {
h = tos;
2008-06-07 05:27:43 +02:00
p0 = cdr(A);
while (iscons(p0)) {
if (!equal(car(p0), F)) {
push(car(p0));
eval(); // force expansion of (x+1)^2, f.e.
}
p0 = cdr(p0);
2008-06-07 03:15:47 +02:00
}
multiply_all(tos - h);
} else
push_integer(1);
T = pop();
}
2008-04-14 02:20:00 +02:00
2008-06-07 04:54:52 +02:00
// Returns the expansion coefficient vector B.
2008-04-13 23:36:03 +02:00
2008-06-07 03:15:47 +02:00
void
2008-06-07 18:07:10 +02:00
expand_get_B(void)
2008-06-07 03:15:47 +02:00
{
2008-06-07 04:54:52 +02:00
int i, n;
if (!istensor(C))
return;
n = C->u.tensor->dim[0];
2008-06-07 03:15:47 +02:00
T = alloc_tensor(n);
T->u.tensor->ndim = 1;
T->u.tensor->dim[0] = n;
for (i = 0; i < n; i++) {
push(B);
2008-04-13 23:36:03 +02:00
push(X);
2008-06-07 03:15:47 +02:00
push_integer(i);
2008-04-13 23:36:03 +02:00
power();
2008-06-07 03:15:47 +02:00
divide();
push(X);
filter();
T->u.tensor->elem[i] = pop();
2008-04-13 23:36:03 +02:00
}
2008-06-07 04:54:52 +02:00
B = T;
2008-04-12 05:28:17 +02:00
}
2008-06-07 04:54:52 +02:00
// Returns the expansion fractions in A.
2008-04-12 05:28:17 +02:00
void
2008-06-07 18:07:10 +02:00
expand_get_A(void)
2008-04-12 05:28:17 +02:00
{
2008-06-07 04:54:52 +02:00
int h, i, n;
if (!istensor(C)) {
push(A);
reciprocate();
A = pop();
return;
}
h = tos;
2008-06-07 03:15:47 +02:00
if (car(A) == symbol(MULTIPLY)) {
T = cdr(A);
while (iscons(T)) {
F = car(T);
2008-06-07 18:07:10 +02:00
expand_get_AF();
2008-06-07 03:15:47 +02:00
T = cdr(T);
}
} else {
F = A;
2008-06-07 18:07:10 +02:00
expand_get_AF();
2008-06-07 03:15:47 +02:00
}
n = tos - h;
T = alloc_tensor(n);
T->u.tensor->ndim = 1;
T->u.tensor->dim[0] = n;
for (i = 0; i < n; i++)
T->u.tensor->elem[i] = stack[h + i];
tos = h;
2008-06-07 04:54:52 +02:00
A = T;
2008-04-12 05:28:17 +02:00
}
2008-04-13 23:36:03 +02:00
2008-04-14 02:20:00 +02:00
void
2008-06-07 18:07:10 +02:00
expand_get_AF(void)
2008-06-07 03:15:47 +02:00
{ int d, i, j, n = 1;
if (!find(F, X))
return;
2008-06-07 18:07:10 +02:00
if (car(F) == symbol(POWER)) {
2008-06-07 03:15:47 +02:00
push(caddr(F));
n = pop_integer();
F = cadr(F);
}
push(F);
2008-04-14 02:20:00 +02:00
push(X);
2008-06-07 03:15:47 +02:00
degree();
d = pop_integer();
for (i = n; i > 0; i--) {
for (j = 0; j < d; j++) {
push(F);
push_integer(i);
power();
reciprocate();
push(X);
push_integer(j);
power();
multiply();
}
}
2008-04-14 02:20:00 +02:00
}
2008-04-13 23:36:03 +02:00
#if SELFTEST
static char *s[] = {
2008-04-20 00:36:37 +02:00
// general cases
2008-04-13 23:36:03 +02:00
"expand(1/(x+1)/(x+2))",
"1/(x+1)-1/(x+2)",
"expand((2x^3-x+2)/(x^2-2x+1))",
"4+2*x+5/(x-1)+3/(x^2-2*x+1)",
2008-04-14 02:20:00 +02:00
"expand(1/x^2/(x-1))",
"-1/(x^2)-1/x+1/(x-1)",
2008-04-20 00:36:37 +02:00
"p=5s+2",
"",
"q=(s+1)(s+2)^2",
"",
"expand(p/q)",
"-3/(s+1)+3/(s+2)+8/(s^2+4*s+4)",
2008-04-13 23:36:03 +02:00
// ensure denominators are expanded (result seems preferable that way)
"q=(x-1)(x-2)^3",
"",
"expand(1/q)",
"1/(x^3-6*x^2+12*x-8)+1/(x-2)-1/(x-1)-1/(x^2-4*x+4)",
2008-04-14 04:35:01 +02:00
// fractional poles
"expand(1/(x+1/2)/(x+1/3))",
2008-06-07 05:27:43 +02:00
"-12/(2*x+1)+18/(3*x+1)",
2008-04-19 03:28:05 +02:00
// expand tensor
"f=1/(x+1)/(x+2)",
"",
"g=1/(x+1)-1/(x+2)",
"",
"expand(((f,f),(f,f)))-((g,g),(g,g))",
"((0,0),(0,0))",
2008-04-19 04:24:34 +02:00
// denominator normalized?
"expand(1/(1+1/x))",
"1-1/(x+1)",
2008-04-19 05:05:28 +02:00
// poles at zero
"expand(1/x/(x+1))",
"1/x-1/(x+1)",
"expand(1/x^2/(x+1))",
"x^(-2)-1/x+1/(x+1)",
// other corner cases
"expand(1/x)",
"1/x",
"expand(1/x^2)",
"x^(-2)",
2008-05-03 20:18:06 +02:00
"expand(1/(x^2-4x+4))",
"1/(x^2-4*x+4)",
2008-04-13 23:36:03 +02:00
};
void
test_expand(void)
{
test(__FILE__, s, sizeof s / sizeof (char *));
}
#endif