add.cpp

   1 /* Symbolic addition
   2
   3         Terms in a sum are combined if they are identical modulo rational
   4         coefficients.
   5
   6         For example, A + 2A becomes 3A.
   7
   8         However, the sum A + sqrt(2) A is not modified.
   9
  10         Combining terms can lead to second-order effects.
  11
  12         For example, consider the case of
  13
  14                 1/sqrt(2) A + 3/sqrt(2) A + sqrt(2) A
  15
  16         The first two terms are combined to yield 2 sqrt(2) A.
  17
  18         This result can now be combined with the third term to yield
  19
  20                 3 sqrt(2) A
  21 */
  22
  23 #include "stdafx.h"
  24 #include "defs.h"
  25
  26 static int flag;
  27
  28 void
  29 eval_add(void)
  30 {
  31         int h = tos;
  32         p1 = cdr(p1);
  33         while (iscons(p1)) {
  34                 push(car(p1));
  35                 eval();
  36                 p2 = pop();
  37                 push_terms(p2);
  38                 p1 = cdr(p1);
  39         }
  40         add_terms(tos - h);
  41 }
  42
  43 /* Add n terms, returns one expression on the stack. */
  44
  45 void
  46 add_terms(int n)
  47 {
  48         int i, h;
  49         U **s;
  50
  51         h = tos - n;
  52
  53         s = stack + h;
  54
  55         /* ensure no infinite loop, use "for" */
  56
  57         for (i = 0; i < 10; i++) {
  58
  59                 if (n < 2)
  60                         break;
  61
  62                 flag = 0;
  63
  64                 qsort(s, n, sizeof (U *), cmp_terms);
  65
  66                 if (flag == 0)
  67                         break;
  68
  69                 n = combine_terms(s, n);
  70         }
  71
  72         tos = h + n;
  73
  74         switch (n) {
  75         case 0:
  76                 push_integer(0);
  77                 break;
  78         case 1:
  79                 break;
  80         default:
  81                 list(n);
  82                 p1 = pop();
  83                 push_symbol(ADD);
  84                 push(p1);
  85                 cons();
  86                 break;
  87         }
  88 }
  89
  90 /* Compare terms for order, clobbers p1 and p2. */
  91
  92 int
  93 cmp_terms(const void *q1, const void *q2)
  94 {
  95         int i, t;
  96
  97         p1 = *((U **) q1);
  98         p2 = *((U **) q2);
  99
 100         /* numbers can be combined */
 101
 102         if (isnum(p1) && isnum(p2)) {
 103                 flag = 1;
 104                 return 0;
 105         }
 106
 107         /* congruent tensors can be combined */
 108
 109         if (istensor(p1) && istensor(p2)) {
 110                 if (p1->u.tensor->ndim < p2->u.tensor->ndim)
 111                         return -1;
 112                 if (p1->u.tensor->ndim > p2->u.tensor->ndim)
 113                         return 1;
 114                 for (i = 0; i < p1->u.tensor->ndim; i++) {
 115                         if (p1->u.tensor->dim[i] < p2->u.tensor->dim[i])
 116                                 return -1;
 117                         if (p1->u.tensor->dim[i] > p2->u.tensor->dim[i])
 118                                 return 1;
 119                 }
 120                 flag = 1;
 121                 return 0;
 122         }
 123
 124         if (car(p1) == symbol(MULTIPLY)) {
 125                 p1 = cdr(p1);
 126                 if (isnum(car(p1))) {
 127                         p1 = cdr(p1);
 128                         if (cdr(p1) == symbol(NIL))
 129                                 p1 = car(p1);
 130                 }
 131         }
 132
 133         if (car(p2) == symbol(MULTIPLY)) {
 134                 p2 = cdr(p2);
 135                 if (isnum(car(p2))) {
 136                         p2 = cdr(p2);
 137                         if (cdr(p2) == symbol(NIL))
 138                                 p2 = car(p2);
 139                 }
 140         }
 141
 142         t = cmp_expr(p1, p2);
 143
 144         if (t == 0)
 145                 flag = 1;
 146
 147         return t;
 148 }
 149
 150 /* Compare adjacent terms in s[] and combine if possible.
 151
 152         Returns the number of terms remaining in s[].
 153
 154         n       number of terms in s[] initially
 155 */
 156
 157 int
 158 combine_terms(U **s, int n)
 159 {
 160         int i, j, t;
 161
 162         for (i = 0; i < n - 1; i++) {
 163                 check_esc_flag();
 164
 165                 p3 = s[i];
 166                 p4 = s[i + 1];
 167
 168                 if (istensor(p3) && istensor(p4)) {
 169                         push(p3);
 170                         push(p4);
 171                         tensor_plus_tensor();
 172                         p1 = pop();
 173                         if (p1 != symbol(NIL)) {
 174                                 s[i] = p1;
 175                                 for (j = i + 1; j < n - 1; j++)
 176                                         s[j] = s[j + 1];
 177                                 n--;
 178                                 i--;
 179                         }
 180                         continue;
 181                 }
 182
 183                 if (istensor(p3) || istensor(p4))
 184                         continue;
 185
 186                 if (isnum(p3) && isnum(p4)) {
 187                         push(p3);
 188                         push(p4);
 189                         add_numbers();
 190                         p1 = pop();
 191                         if (iszero(p1)) {
 192                                 for (j = i; j < n - 2; j++)
 193                                         s[j] = s[j + 2];
 194                                 n -= 2;
 195                         } else {
 196                                 s[i] = p1;
 197                                 for (j = i + 1; j < n - 1; j++)
 198                                         s[j] = s[j + 1];
 199                                 n--;
 200                         }
 201                         i--;
 202                         continue;
 203                 }
 204
 205                 if (isnum(p3) || isnum(p4))
 206                         continue;
 207
 208                 p1 = one;
 209                 p2 = one;
 210
 211                 t = 0;
 212
 213                 if (car(p3) == symbol(MULTIPLY)) {
 214                         p3 = cdr(p3);
 215                         t = 1; /* p3 is now denormal */
 216                         if (isnum(car(p3))) {
 217                                 p1 = car(p3);
 218                                 p3 = cdr(p3);
 219                                 if (cdr(p3) == symbol(NIL)) {
 220                                         p3 = car(p3);
 221                                         t = 0;
 222                                 }
 223                         }
 224                 }
 225
 226                 if (car(p4) == symbol(MULTIPLY)) {
 227                         p4 = cdr(p4);
 228                         if (isnum(car(p4))) {
 229                                 p2 = car(p4);
 230                                 p4 = cdr(p4);
 231                                 if (cdr(p4) == symbol(NIL))
 232                                         p4 = car(p4);
 233                         }
 234                 }
 235
 236                 if (!equal(p3, p4))
 237                         continue;
 238
 239                 push(p1);
 240                 push(p2);
 241                 add_numbers();
 242
 243                 p1 = pop();
 244
 245                 if (iszero(p1)) {
 246                         for (j = i; j < n - 2; j++)
 247                                 s[j] = s[j + 2];
 248                         n -= 2;
 249                         i--;
 250                         continue;
 251                 }
 252
 253                 push(p1);
 254
 255                 if (t) {
 256                         push(symbol(MULTIPLY));
 257                         push(p3);
 258                         cons();
 259                 } else
 260                         push(p3);
 261
 262                 multiply();
 263
 264                 s[i] = pop();
 265
 266                 for (j = i + 1; j < n - 1; j++)
 267                         s[j] = s[j + 1];
 268
 269                 n--;
 270                 i--;
 271         }
 272
 273         return n;
 274 }
 275
 276 void
 277 push_terms(U *p)
 278 {
 279         if (car(p) == symbol(ADD)) {
 280                 p = cdr(p);
 281                 while (iscons(p)) {
 282                         push(car(p));
 283                         p = cdr(p);
 284                 }
 285         } else if (!iszero(p))
 286                 push(p);
 287 }
 288
 289 /* add two expressions */
 290
 291 void
 292 add()
 293 {
 294         int h;
 295         save();
 296         p2 = pop();
 297         p1 = pop();
 298         h = tos;
 299         push_terms(p1);
 300         push_terms(p2);
 301         add_terms(tos - h);
 302         restore();
 303 }
 304
 305 void
 306 add_all(int k)
 307 {
 308         int h, i;
 309         U **s;
 310         save();
 311         s = stack + tos - k;
 312         h = tos;
 313         for (i = 0; i < k; i++)
 314                 push_terms(s[i]);
 315         add_terms(tos - h);
 316         p1 = pop();
 317         tos -= k;
 318         push(p1);
 319         restore();
 320 }
 321
 322 void
 323 subtract(void)
 324 {
 325         negate();
 326         add();
 327 }