rationalize.cpp

   1 #include "stdafx.h"
   2 #include "defs.h"
   3 #define DEBUG 0
   4 static void __rationalize_tensor(void);
   5 static void multiply_denominators(U *);
   6 static void multiply_denominators_term(U *);
   7 static void multiply_denominators_factor(U *);
   8 static void __lcm(void);
   9
  10 void
  11 eval_rationalize(void)
  12 {
  13         push(cadr(p1));
  14         eval();
  15         rationalize();
  16 }
  17
  18 void
  19 rationalize(void)
  20 {
  21         int x = expanding;
  22         save();
  23         yyrationalize();
  24         restore();
  25         expanding = x;
  26 }
  27
  28 void
  29 yyrationalize(void)
  30 {
  31         p1 = pop();
  32
  33         if (istensor(p1)) {
  34                 __rationalize_tensor();
  35                 return;
  36         }
  37
  38         expanding = 0;
  39
  40         if (car(p1) != symbol(ADD)) {
  41                 push(p1);
  42                 return;
  43         }
  44
  45         // get common denominator
  46
  47         push(one);
  48         multiply_denominators(p1);
  49         p2 = pop();
  50
  51         // multiply each term by common denominator
  52
  53         push(zero);
  54         p3 = cdr(p1);
  55         while (iscons(p3)) {
  56                 push(p2);
  57                 push(car(p3));
  58                 multiply();
  59                 add();
  60                 p3 = cdr(p3);
  61         }
  62
  63         // collect common factors
  64
  65         Condense();
  66
  67         // divide by common denominator
  68
  69         push(p2);
  70         divide();
  71 }
  72
  73 static void
  74 multiply_denominators(U *p)
  75 {
  76         if (car(p) == symbol(ADD)) {
  77                 p = cdr(p);
  78                 while (iscons(p)) {
  79                         multiply_denominators_term(car(p));
  80                         p = cdr(p);
  81                 }
  82         } else
  83                 multiply_denominators_term(p);
  84 }
  85
  86 static void
  87 multiply_denominators_term(U *p)
  88 {
  89         if (car(p) == symbol(MULTIPLY)) {
  90                 p = cdr(p);
  91                 while (iscons(p)) {
  92                         multiply_denominators_factor(car(p));
  93                         p = cdr(p);
  94                 }
  95         } else
  96                 multiply_denominators_factor(p);
  97 }
  98
  99 static void
 100 multiply_denominators_factor(U *p)
 101 {
 102         if (car(p) != symbol(POWER))
 103                 return;
 104
 105         push(p);
 106
 107         p = caddr(p);
 108
 109         // like x^(-2) ?
 110
 111         if (isnegativenumber(p)) {
 112                 inverse();
 113                 __lcm();
 114                 return;
 115         }
 116
 117         // like x^(-a) ?
 118
 119         if (car(p) == symbol(MULTIPLY) && isnegativenumber(cadr(p))) {
 120                 inverse();
 121                 __lcm();
 122                 return;
 123         }
 124
 125         // no match
 126
 127         pop();
 128 }
 129
 130 static void
 131 __rationalize_tensor(void)
 132 {
 133         int i, n;
 134
 135         push(p1);
 136
 137         eval(); // makes a copy
 138
 139         p1 = pop();
 140
 141         if (!istensor(p1)) { // might be zero
 142                 push(p1);
 143                 return;
 144         }
 145
 146         n = p1->u.tensor->nelem;
 147
 148         for (i = 0; i < n; i++) {
 149                 push(p1->u.tensor->elem[i]);
 150                 rationalize();
 151                 p1->u.tensor->elem[i] = pop();
 152         }
 153
 154         push(p1);
 155 }
 156
 157 static void
 158 __lcm(void)
 159 {
 160         save();
 161
 162         p1 = pop();
 163         p2 = pop();
 164
 165         push(p1);
 166         push(p2);
 167         multiply();
 168         push(p1);
 169         push(p2);
 170         gcd();
 171         divide();
 172
 173         restore();
 174 }