simfac.cpp

   1 /* Simplify factorials
   2
   3 The following script
   4
   5         F(n,k) = k binomial(n,k)
   6         (F(n,k) + F(n,k-1)) / F(n+1,k)
   7
   8 generates
   9
  10        k! n!             n! (1 - k + n)!              k! n!
  11  -------------------- + -------------------- - ----------------------
  12  (-1 + k)! (1 + n)!     (1 + n)! (-k + n)!     k (-1 + k)! (1 + n)!
  13
  14 Simplify each term to get
  15
  16    k       1 - k + n       1
  17 ------- + ----------- - -------
  18  1 + n       1 + n       1 + n
  19
  20 Then simplify the sum to get
  21
  22    n
  23 -------
  24  1 + n
  25
  26 */
  27
  28 #include "stdafx.h"
  29 #include "defs.h"
  30 void simfac(void);
  31 static void simfac_term(void);
  32 static int yysimfac(int);
  33
  34 // simplify factorials term-by-term
  35
  36 void
  37 eval_simfac(void)
  38 {
  39         push(cadr(p1));
  40         eval();
  41         simfac();
  42 }
  43
  44 #if 1
  45
  46 void
  47 simfac(void)
  48 {
  49         int h;
  50         save();
  51         p1 = pop();
  52         if (car(p1) == symbol(ADD)) {
  53                 h = tos;
  54                 p1 = cdr(p1);
  55                 while (p1 != symbol(NIL)) {
  56                         push(car(p1));
  57                         simfac_term();
  58                         p1 = cdr(p1);
  59                 }
  60                 add_all(tos - h);
  61         } else {
  62                 push(p1);
  63                 simfac_term();
  64         }
  65         restore();
  66 }
  67
  68 #else
  69
  70 void
  71 simfac(void)
  72 {
  73         int h;
  74         save();
  75         p1 = pop();
  76         if (car(p1) == symbol(ADD)) {
  77                 h = tos;
  78                 p1 = cdr(p1);
  79                 while (p1 != symbol(NIL)) {
  80                         push(car(p1));
  81                         simfac_term();
  82                         p1 = cdr(p1);
  83                 }
  84                 addk(tos - h);
  85                 p1 = pop();
  86                 if (find(p1, symbol(FACTORIAL))) {
  87                         push(p1);
  88                         if (car(p1) == symbol(ADD)) {
  89                                 Condense();
  90                                 simfac_term();
  91                         }
  92                 }
  93         } else {
  94                 push(p1);
  95                 simfac_term();
  96         }
  97         restore();
  98 }
  99
 100 #endif
 101
 102 static void
 103 simfac_term(void)
 104 {
 105         int h;
 106
 107         save();
 108
 109         p1 = pop();
 110
 111         // if not a product of factors then done
 112
 113         if (car(p1) != symbol(MULTIPLY)) {
 114                 push(p1);
 115                 restore();
 116                 return;
 117         }
 118
 119         // push all factors
 120
 121         h = tos;
 122         p1 = cdr(p1);
 123         while (p1 != symbol(NIL)) {
 124                 push(car(p1));
 125                 p1 = cdr(p1);
 126         }
 127
 128         // keep trying until no more to do
 129
 130         while (yysimfac(h))
 131                 ;
 132
 133         multiply_all_noexpand(tos - h);
 134
 135         restore();
 136 }
 137
 138 // try all pairs of factors
 139
 140 static int
 141 yysimfac(int h)
 142 {
 143         int i, j;
 144
 145         for (i = h; i < tos; i++) {
 146                 p1 = stack[i];
 147                 for (j = h; j < tos; j++) {
 148                         if (i == j)
 149                                 continue;
 150                         p2 = stack[j];
 151
 152                         //      n! / n          ->      (n - 1)!
 153
 154                         if (car(p1) == symbol(FACTORIAL)
 155                         && car(p2) == symbol(POWER)
 156                         && isminusone(caddr(p2))
 157                         && equal(cadr(p1), cadr(p2))) {
 158                                 push(cadr(p1));
 159                                 push(one);
 160                                 subtract();
 161                                 factorial();
 162                                 stack[i] = pop();
 163                                 stack[j] = one;
 164                                 return 1;
 165                         }
 166
 167                         //      n / n!          ->      1 / (n - 1)!
 168
 169                         if (car(p2) == symbol(POWER)
 170                         && isminusone(caddr(p2))
 171                         && caadr(p2) == symbol(FACTORIAL)
 172                         && equal(p1, cadadr(p2))) {
 173                                 push(p1);
 174                                 push_integer(-1);
 175                                 add();
 176                                 factorial();
 177                                 reciprocate();
 178                                 stack[i] = pop();
 179                                 stack[j] = one;
 180                                 return 1;
 181                         }
 182
 183                         //      (n + 1) n!      ->      (n + 1)!
 184
 185                         if (car(p2) == symbol(FACTORIAL)) {
 186                                 push(p1);
 187                                 push(cadr(p2));
 188                                 subtract();
 189                                 p3 = pop();
 190                                 if (isplusone(p3)) {
 191                                         push(p1);
 192                                         factorial();
 193                                         stack[i] = pop();
 194                                         stack[j] = one;
 195                                         return 1;
 196                                 }
 197                         }
 198
 199                         //      1 / ((n + 1) n!)        ->      1 / (n + 1)!
 200
 201                         if (car(p1) == symbol(POWER)
 202                         && isminusone(caddr(p1))
 203                         && car(p2) == symbol(POWER)
 204                         && isminusone(caddr(p2))
 205                         && caadr(p2) == symbol(FACTORIAL)) {
 206                                 push(cadr(p1));
 207                                 push(cadr(cadr(p2)));
 208                                 subtract();
 209                                 p3 = pop();
 210                                 if (isplusone(p3)) {
 211                                         push(cadr(p1));
 212                                         factorial();
 213                                         reciprocate();
 214                                         stack[i] = pop();
 215                                         stack[j] = one;
 216                                         return 1;
 217                                 }
 218                         }
 219
 220                         //      (n + 1)! / n!   ->      n + 1
 221
 222                         //      n! / (n + 1)!   ->      1 / (n + 1)
 223
 224                         if (car(p1) == symbol(FACTORIAL)
 225                         && car(p2) == symbol(POWER)
 226                         && isminusone(caddr(p2))
 227                         && caadr(p2) == symbol(FACTORIAL)) {
 228                                 push(cadr(p1));
 229                                 push(cadr(cadr(p2)));
 230                                 subtract();
 231                                 p3 = pop();
 232                                 if (isplusone(p3)) {
 233                                         stack[i] = cadr(p1);
 234                                         stack[j] = one;
 235                                         return 1;
 236                                 }
 237                                 if (isminusone(p3)) {
 238                                         push(cadr(cadr(p2)));
 239                                         reciprocate();
 240                                         stack[i] = pop();
 241                                         stack[j] = one;
 242                                         return 1;
 243                                 }
 244                                 if (equaln(p3, 2)) {
 245                                         stack[i] = cadr(p1);
 246                                         push(cadr(p1));
 247                                         push_integer(-1);
 248                                         add();
 249                                         stack[j] = pop();
 250                                         return 1;
 251                                 }
 252                                 if (equaln(p3, -2)) {
 253                                         push(cadr(cadr(p2)));
 254                                         reciprocate();
 255                                         stack[i] = pop();
 256                                         push(cadr(cadr(p2)));
 257                                         push_integer(-1);
 258                                         add();
 259                                         reciprocate();
 260                                         stack[j] = pop();
 261                                         return 1;
 262                                 }
 263                         }
 264                 }
 265         }
 266
 267         return 0;
 268 }