Eigenmath/rationalize.cpp

#include "stdafx.h"
#include "defs.h"
#define DEBUG 0
static void __rationalize_tensor(void);
static void multiply_denominators(U *);
static void multiply_denominators_term(U *);
static void multiply_denominators_factor(U *);
static void __lcm(void);

void
eval_rationalize(void)
{
	push(cadr(p1));
	eval();
	rationalize();
}

void
rationalize(void)
{
	int x = expanding;
	save();
	yyrationalize();
	restore();
	expanding = x;
}

void
yyrationalize(void)
{
	p1 = pop();

	if (istensor(p1)) {
		__rationalize_tensor();
		return;
	}

	expanding = 0;

	if (car(p1) != symbol(ADD)) {
		push(p1);
		return;
	}

	// get common denominator

	push(one);
	multiply_denominators(p1);
	p2 = pop();

	// multiply each term by common denominator

	push(zero);
	p3 = cdr(p1);
	while (iscons(p3)) {
		push(p2);
		push(car(p3));
		multiply();
		add();
		p3 = cdr(p3);
	}

	// collect common factors

	Condense();

	// divide by common denominator

	push(p2);
	divide();
}

static void
multiply_denominators(U *p)
{
	if (car(p) == symbol(ADD)) {
		p = cdr(p);
		while (iscons(p)) {
			multiply_denominators_term(car(p));
			p = cdr(p);
		}
	} else
		multiply_denominators_term(p);
}

static void
multiply_denominators_term(U *p)
{
	if (car(p) == symbol(MULTIPLY)) {
		p = cdr(p);
		while (iscons(p)) {
			multiply_denominators_factor(car(p));
			p = cdr(p);
		}
	} else
		multiply_denominators_factor(p);
}

static void
multiply_denominators_factor(U *p)
{
	if (car(p) != symbol(POWER))
		return;

	push(p);

	p = caddr(p);

	// like x^(-2) ?

	if (isnegativenumber(p)) {
		inverse();
		__lcm();
		return;
	}

	// like x^(-a) ?

	if (car(p) == symbol(MULTIPLY) && isnegativenumber(cadr(p))) {
		inverse();
		__lcm();
		return;
	}

	// no match

	pop();
}

static void
__rationalize_tensor(void)
{
	int i, n;

	push(p1);

	eval(); // makes a copy

	p1 = pop();

	if (!istensor(p1)) { // might be zero
		push(p1);
		return;
	}

	n = p1->u.tensor->nelem;

	for (i = 0; i < n; i++) {
		push(p1->u.tensor->elem[i]);
		rationalize();
		p1->u.tensor->elem[i] = pop();
	}

	push(p1);
}

static void
__lcm(void)
{
	save();

	p1 = pop();
	p2 = pop();

	push(p1);
	push(p2);
	multiply();
	push(p1);
	push(p2);
	gcd();
	divide();

	restore();
}