123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175 |
- #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();
- }
|