simplify.cpp

   1 #include "stdafx.h"
   2 #include "defs.h"
   3
   4 extern int trigmode;
   5 static void simplify_tensor(void);
   6 static int count(U *);
   7 static int nterms(U *);
   8 static void f1(void);
   9 static void f2(void);
  10 static void f3(void);
  11 static void f4(void);
  12 static void f5(void);
  13 static void f9(void);
  14
  15 void
  16 eval_simplify(void)
  17 {
  18         push(cadr(p1));
  19         eval();
  20         simplify();
  21 }
  22
  23 void
  24 simplify(void)
  25 {
  26         save();
  27         simplify_main();
  28         restore();
  29 }
  30
  31 void
  32 simplify_main(void)
  33 {
  34         p1 = pop();
  35
  36         if (istensor(p1)) {
  37                 simplify_tensor();
  38                 return;
  39         }
  40
  41         if (find(p1, symbol(FACTORIAL))) {
  42                 push(p1);
  43                 simfac();
  44                 p2 = pop();
  45                 push(p1);
  46                 rationalize();
  47                 simfac();
  48                 p3 = pop();
  49                 if (count(p2) < count(p3))
  50                         p1 = p2;
  51                 else
  52                         p1 = p3;
  53         }
  54
  55         f1();
  56         f2();
  57         f3();
  58         f4();
  59         f5();
  60         f9();
  61
  62         push(p1);
  63 }
  64
  65 static void
  66 simplify_tensor(void)
  67 {
  68         int i;
  69         p2 = alloc_tensor(p1->u.tensor->nelem);
  70         p2->u.tensor->ndim = p1->u.tensor->ndim;
  71         for (i = 0; i < p1->u.tensor->ndim; i++)
  72                 p2->u.tensor->dim[i] = p1->u.tensor->dim[i];
  73         for (i = 0; i < p1->u.tensor->nelem; i++) {
  74                 push(p1->u.tensor->elem[i]);
  75                 simplify();
  76                 p2->u.tensor->elem[i] = pop();
  77         }
  78         if (iszero(p2))
  79                 p2 = zero; // null tensor becomes scalar zero
  80         push(p2);
  81 }
  82
  83 static int
  84 count(U *p)
  85 {
  86         int n;
  87         if (iscons(p)) {
  88                 n = 0;
  89                 while (iscons(p)) {
  90                         n += count(car(p)) + 1;
  91                         p = cdr(p);
  92                 }
  93         } else
  94                 n = 1;
  95         return n;
  96 }
  97
  98 // try rationalizing
  99
 100 static void
 101 f1(void)
 102 {
 103         if (car(p1) != symbol(ADD))
 104                 return;
 105         push(p1);
 106         rationalize();
 107         p2 = pop();
 108         if (count(p2) < count(p1))
 109                 p1 = p2;
 110 }
 111
 112 // try condensing
 113
 114 static void
 115 f2(void)
 116 {
 117         if (car(p1) != symbol(ADD))
 118                 return;
 119         push(p1);
 120         Condense();
 121         p2 = pop();
 122         if (count(p2) <= count(p1))
 123                 p1 = p2;
 124 }
 125
 126 // this simplifies forms like (A-B) / (B-A)
 127
 128 static void
 129 f3(void)
 130 {
 131         push(p1);
 132         rationalize();
 133         negate();
 134         rationalize();
 135         negate();
 136         rationalize();
 137         p2 = pop();
 138         if (count(p2) < count(p1))
 139                 p1 = p2;
 140 }
 141
 142 // try expanding denominators
 143
 144 static void
 145 f4(void)
 146 {
 147         if (iszero(p1))
 148                 return;
 149         push(p1);
 150         rationalize();
 151         inverse();
 152         rationalize();
 153         inverse();
 154         rationalize();
 155         p2 = pop();
 156         if (count(p2) < count(p1))
 157                 p1 = p2;
 158 }
 159
 160 // simplifies trig forms
 161
 162 void
 163 simplify_trig(void)
 164 {
 165         save();
 166         p1 = pop();
 167         f5();
 168         push(p1);
 169         restore();
 170 }
 171
 172 static void
 173 f5(void)
 174 {
 175         if (find(p1, symbol(SIN)) == 0 && find(p1, symbol(COS)) == 0)
 176                 return;
 177
 178         p2 = p1;
 179
 180         trigmode = 1;
 181         push(p2);
 182         eval();
 183         p3 = pop();
 184
 185         trigmode = 2;
 186         push(p2);
 187         eval();
 188         p4 = pop();
 189
 190         trigmode = 0;
 191
 192         if (count(p4) < count(p3) || nterms(p4) < nterms(p3))
 193                 p3 = p4;
 194
 195         if (count(p3) < count(p1) || nterms(p3) < nterms(p1))
 196                 p1 = p3;
 197 }
 198
 199 // if it's a sum then try to simplify each term
 200
 201 static void
 202 f9(void)
 203 {
 204         if (car(p1) != symbol(ADD))
 205                 return;
 206         push_integer(0);
 207         p2 = cdr(p1);
 208         while (iscons(p2)) {
 209                 push(car(p2));
 210                 simplify();
 211                 add();
 212                 p2 = cdr(p2);
 213         }
 214         p2 = pop();
 215         if (count(p2) < count(p1))
 216                 p1 = p2;
 217 }
 218
 219 static int
 220 nterms(U *p)
 221 {
 222         if (car(p) != symbol(ADD))
 223                 return 1;
 224         else
 225                 return length(p) - 1;
 226 }
 227
 228 #if SELFTEST
 229
 230 static char *s[] = {
 231
 232         "simplify(A)",
 233         "A",
 234
 235         "simplify(A+B)",
 236         "A+B",
 237
 238         "simplify(A B)",
 239         "A*B",
 240
 241         "simplify(A^B)",
 242         "A^B",
 243
 244         "simplify(A/(A+B)+B/(A+B))",
 245         "1",
 246
 247         "simplify((A-B)/(B-A))",
 248         "-1",
 249
 250         "A=((A11,A12),(A21,A22))",
 251         "",
 252
 253         "simplify(det(A) inv(A) - adj(A))",
 254         "0",
 255
 256         "A=quote(A)",
 257         "",
 258
 259         // this shows need for <= in try_factoring
 260
 261 //      "x*(1+a)",
 262 //      "x+a*x",
 263
 264 //      "simplify(last)",
 265 //      "x*(1+a)",
 266
 267         "simplify(-3 exp(-1/3 r + i phi) cos(theta) / sin(theta)"
 268         " + 3 exp(-1/3 r + i phi) cos(theta) sin(theta)"
 269         " + 3 exp(-1/3 r + i phi) cos(theta)^3 / sin(theta))",
 270         "0",
 271
 272         "simplify((A^2 C^2 + A^2 D^2 + B^2 C^2 + B^2 D^2)/(A^2+B^2)/(C^2+D^2))",
 273         "1",
 274
 275         "simplify(d(arctan(y/x),y))",
 276         "x/(x^2+y^2)",
 277
 278         "simplify(d(arctan(y/x),x))",
 279         "-y/(x^2+y^2)",
 280
 281         "simplify(1-sin(x)^2)",
 282         "cos(x)^2",
 283
 284         "simplify(1-cos(x)^2)",
 285         "sin(x)^2",
 286
 287         "simplify(sin(x)^2-1)",
 288         "-cos(x)^2",
 289
 290         "simplify(cos(x)^2-1)",
 291         "-sin(x)^2",
 292 /*
 293         "simfac(n!/n)-(n-1)!",
 294         "0",
 295
 296         "simfac(n/n!)-1/(n-1)!",
 297         "0",
 298
 299         "simfac(rationalize((n+k+1)/(n+k+1)!))-1/(n+k)!",
 300         "0",
 301
 302         "simfac(condense((n+1)*n!))-(n+1)!",
 303         "0",
 304
 305         "simfac(1/((n+1)*n!))-1/(n+1)!",
 306         "0",
 307
 308         "simfac((n+1)!/n!)-n-1",
 309         "0",
 310
 311         "simfac(n!/(n+1)!)-1/(n+1)",
 312         "0",
 313
 314         "simfac(binomial(n+1,k)/binomial(n,k))",
 315         "(1+n)/(1-k+n)",
 316
 317         "simfac(binomial(n,k)/binomial(n+1,k))",
 318         "(1-k+n)/(1+n)",
 319
 320         "F(nn,kk)=kk*binomial(nn,kk)",
 321         "",
 322
 323         "simplify(simfac((F(n,k)+F(n,k-1))/F(n+1,k))-n/(n+1))",
 324         "0",
 325
 326         "F=quote(F)",
 327         "",
 328 */
 329         "simplify(n!/n)-(n-1)!",
 330         "0",
 331
 332         "simplify(n/n!)-1/(n-1)!",
 333         "0",
 334
 335         "simplify(rationalize((n+k+1)/(n+k+1)!))-1/(n+k)!",
 336         "0",
 337
 338         "simplify(condense((n+1)*n!))-(n+1)!",
 339         "0",
 340
 341         "simplify(1/((n+1)*n!))-1/(n+1)!",
 342         "0",
 343
 344         "simplify((n+1)!/n!)-n-1",
 345         "0",
 346
 347         "simplify(n!/(n+1)!)-1/(n+1)",
 348         "0",
 349
 350         "simplify(binomial(n+1,k)/binomial(n,k))",
 351         "(1+n)/(1-k+n)",
 352
 353         "simplify(binomial(n,k)/binomial(n+1,k))",
 354         "(1-k+n)/(1+n)",
 355
 356         "F(nn,kk)=kk*binomial(nn,kk)",
 357         "",
 358
 359         "simplify((F(n,k)+F(n,k-1))/F(n+1,k))-n/(n+1)",
 360         "0",
 361
 362         "F=quote(F)",
 363         "",
 364
 365         "simplify((n+1)/(n+1)!)-1/n!",
 366         "0",
 367
 368         "simplify(a*b+a*c)",
 369         "a*(b+c)",
 370
 371 // Symbol's binding is evaluated, undoing simplify
 372
 373         "x=simplify(a*b+a*c)",
 374         "",
 375
 376         "x",
 377         "a*b+a*c",
 378 };
 379
 380 void
 381 test_simplify(void)
 382 {
 383         test(__FILE__, s, sizeof s / sizeof (char *));
 384 }
 385
 386 #endif