eigenmath/simfac.cpp

269 lines
4.1 KiB
C++

/* Simplify factorials
The following script
F(n,k) = k binomial(n,k)
(F(n,k) + F(n,k-1)) / F(n+1,k)
generates
k! n! n! (1 - k + n)! k! n!
-------------------- + -------------------- - ----------------------
(-1 + k)! (1 + n)! (1 + n)! (-k + n)! k (-1 + k)! (1 + n)!
Simplify each term to get
k 1 - k + n 1
------- + ----------- - -------
1 + n 1 + n 1 + n
Then simplify the sum to get
n
-------
1 + n
*/
#include "stdafx.h"
#include "defs.h"
void simfac(void);
static void simfac_term(void);
static int yysimfac(int);
// simplify factorials term-by-term
void
eval_simfac(void)
{
push(cadr(p1));
eval();
simfac();
}
#if 1
void
simfac(void)
{
int h;
save();
p1 = pop();
if (car(p1) == symbol(ADD)) {
h = tos;
p1 = cdr(p1);
while (p1 != symbol(NIL)) {
push(car(p1));
simfac_term();
p1 = cdr(p1);
}
add_all(tos - h);
} else {
push(p1);
simfac_term();
}
restore();
}
#else
void
simfac(void)
{
int h;
save();
p1 = pop();
if (car(p1) == symbol(ADD)) {
h = tos;
p1 = cdr(p1);
while (p1 != symbol(NIL)) {
push(car(p1));
simfac_term();
p1 = cdr(p1);
}
addk(tos - h);
p1 = pop();
if (find(p1, symbol(FACTORIAL))) {
push(p1);
if (car(p1) == symbol(ADD)) {
Condense();
simfac_term();
}
}
} else {
push(p1);
simfac_term();
}
restore();
}
#endif
static void
simfac_term(void)
{
int h;
save();
p1 = pop();
// if not a product of factors then done
if (car(p1) != symbol(MULTIPLY)) {
push(p1);
restore();
return;
}
// push all factors
h = tos;
p1 = cdr(p1);
while (p1 != symbol(NIL)) {
push(car(p1));
p1 = cdr(p1);
}
// keep trying until no more to do
while (yysimfac(h))
;
multiply_all_noexpand(tos - h);
restore();
}
// try all pairs of factors
static int
yysimfac(int h)
{
int i, j;
for (i = h; i < tos; i++) {
p1 = stack[i];
for (j = h; j < tos; j++) {
if (i == j)
continue;
p2 = stack[j];
// n! / n -> (n - 1)!
if (car(p1) == symbol(FACTORIAL)
&& car(p2) == symbol(POWER)
&& isminusone(caddr(p2))
&& equal(cadr(p1), cadr(p2))) {
push(cadr(p1));
push(one);
subtract();
factorial();
stack[i] = pop();
stack[j] = one;
return 1;
}
// n / n! -> 1 / (n - 1)!
if (car(p2) == symbol(POWER)
&& isminusone(caddr(p2))
&& caadr(p2) == symbol(FACTORIAL)
&& equal(p1, cadadr(p2))) {
push(p1);
push_integer(-1);
add();
factorial();
reciprocate();
stack[i] = pop();
stack[j] = one;
return 1;
}
// (n + 1) n! -> (n + 1)!
if (car(p2) == symbol(FACTORIAL)) {
push(p1);
push(cadr(p2));
subtract();
p3 = pop();
if (isplusone(p3)) {
push(p1);
factorial();
stack[i] = pop();
stack[j] = one;
return 1;
}
}
// 1 / ((n + 1) n!) -> 1 / (n + 1)!
if (car(p1) == symbol(POWER)
&& isminusone(caddr(p1))
&& car(p2) == symbol(POWER)
&& isminusone(caddr(p2))
&& caadr(p2) == symbol(FACTORIAL)) {
push(cadr(p1));
push(cadr(cadr(p2)));
subtract();
p3 = pop();
if (isplusone(p3)) {
push(cadr(p1));
factorial();
reciprocate();
stack[i] = pop();
stack[j] = one;
return 1;
}
}
// (n + 1)! / n! -> n + 1
// n! / (n + 1)! -> 1 / (n + 1)
if (car(p1) == symbol(FACTORIAL)
&& car(p2) == symbol(POWER)
&& isminusone(caddr(p2))
&& caadr(p2) == symbol(FACTORIAL)) {
push(cadr(p1));
push(cadr(cadr(p2)));
subtract();
p3 = pop();
if (isplusone(p3)) {
stack[i] = cadr(p1);
stack[j] = one;
return 1;
}
if (isminusone(p3)) {
push(cadr(cadr(p2)));
reciprocate();
stack[i] = pop();
stack[j] = one;
return 1;
}
if (equaln(p3, 2)) {
stack[i] = cadr(p1);
push(cadr(p1));
push_integer(-1);
add();
stack[j] = pop();
return 1;
}
if (equaln(p3, -2)) {
push(cadr(cadr(p2)));
reciprocate();
stack[i] = pop();
push(cadr(cadr(p2)));
push_integer(-1);
add();
reciprocate();
stack[j] = pop();
return 1;
}
}
}
}
return 0;
}