Eigenmath/simfac.cpp

/* 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;
}