expand.cpp

   1 // Partial fraction expansion
   2 //
   3 // Example
   4 //
   5 //      expand(1/(x^3+x^2),x)
   6 //
   7 //        1      1       1
   8 //      ---- - --- + -------
   9 //        2     x     x + 1
  10 //       x
  11
  12 #include "stdafx.h"
  13 #include "defs.h"
  14
  15 void
  16 eval_expand(void)
  17 {
  18         // 1st arg
  19
  20         push(cadr(p1));
  21         eval();
  22
  23         // 2nd arg
  24
  25         push(caddr(p1));
  26         eval();
  27
  28         p2 = pop();
  29         if (p2 == symbol(NIL))
  30                 guess();
  31         else
  32                 push(p2);
  33
  34         expand();
  35 }
  36
  37 #define A p2
  38 #define B p3
  39 #define C p4
  40 #define F p5
  41 #define P p6
  42 #define Q p7
  43 #define T p8
  44 #define X p9
  45
  46 void
  47 expand(void)
  48 {
  49         save();
  50
  51         X = pop();
  52         F = pop();
  53
  54         if (istensor(F)) {
  55                 expand_tensor();
  56                 restore();
  57                 return;
  58         }
  59
  60         // if sum of terms then sum over the expansion of each term
  61
  62         if (car(F) == symbol(ADD)) {
  63                 push_integer(0);
  64                 p1 = cdr(F);
  65                 while (iscons(p1)) {
  66                         push(car(p1));
  67                         push(X);
  68                         expand();
  69                         add();
  70                         p1 = cdr(p1);
  71                 }
  72                 restore();
  73                 return;
  74         }
  75
  76         // B = numerator
  77
  78         push(F);
  79         numerator();
  80         B = pop();
  81
  82         // A = denominator
  83
  84         push(F);
  85         denominator();
  86         A = pop();
  87
  88         remove_negative_exponents();
  89
  90         // Q = quotient
  91
  92         push(B);
  93         push(A);
  94         push(X);
  95         divpoly();
  96         Q = pop();
  97
  98         // remainder B = B - A * Q
  99
 100         push(B);
 101         push(A);
 102         push(Q);
 103         multiply();
 104         subtract();
 105         B = pop();
 106
 107         // if the remainder is zero then we're done
 108
 109         if (iszero(B)) {
 110                 push(Q);
 111                 restore();
 112                 return;
 113         }
 114
 115         // A = factor(A)
 116
 117         push(A);
 118         push(X);
 119         factorpoly();
 120         A = pop();
 121
 122         expand_get_C();
 123         expand_get_B();
 124         expand_get_A();
 125
 126         if (istensor(C)) {
 127                 push(C);
 128                 inv();
 129                 push(B);
 130                 inner();
 131                 push(A);
 132                 inner();
 133         } else {
 134                 push(B);
 135                 push(C);
 136                 divide();
 137                 push(A);
 138                 multiply();
 139         }
 140
 141         push(Q);
 142         add();
 143
 144         restore();
 145 }
 146
 147 void
 148 expand_tensor(void)
 149 {
 150         int i;
 151         push(F);
 152         copy_tensor();
 153         F = pop();
 154         for (i = 0; i < F->u.tensor->nelem; i++) {
 155                 push(F->u.tensor->elem[i]);
 156                 push(X);
 157                 expand();
 158                 F->u.tensor->elem[i] = pop();
 159         }
 160         push(F);
 161 }
 162
 163 void
 164 remove_negative_exponents(void)
 165 {
 166         int h, i, j, k, n;
 167
 168         h = tos;
 169         factors(A);
 170         factors(B);
 171         n = tos - h;
 172
 173         // find the smallest exponent
 174
 175         j = 0;
 176         for (i = 0; i < n; i++) {
 177                 p1 = stack[h + i];
 178                 if (car(p1) != symbol(POWER))
 179                         continue;
 180                 if (cadr(p1) != X)
 181                         continue;
 182                 push(caddr(p1));
 183                 k = pop_integer();
 184                 if (k == (int) 0x80000000)
 185                         continue;
 186                 if (k < j)
 187                         j = k;
 188         }
 189
 190         tos = h;
 191
 192         if (j == 0)
 193                 return;
 194
 195         // A = A / X^j
 196
 197         push(A);
 198         push(X);
 199         push_integer(-j);
 200         power();
 201         multiply();
 202         A = pop();
 203
 204         // B = B / X^j
 205
 206         push(B);
 207         push(X);
 208         push_integer(-j);
 209         power();
 210         multiply();
 211         B = pop();
 212 }
 213
 214 // Returns the expansion coefficient matrix C.
 215 //
 216 // Example:
 217 //
 218 //       B         1
 219 //      --- = -----------
 220 //       A      2
 221 //             x (x + 1)
 222 //
 223 // We have
 224 //
 225 //       B     Y1     Y2      Y3
 226 //      --- = ---- + ---- + -------
 227 //       A      2     x      x + 1
 228 //             x
 229 //
 230 // Our task is to solve for the unknowns Y1, Y2, and Y3.
 231 //
 232 // Multiplying both sides by A yields
 233 //
 234 //           AY1     AY2      AY3
 235 //      B = ----- + ----- + -------
 236 //            2      x       x + 1
 237 //           x
 238 //
 239 // Let
 240 //
 241 //            A               A                 A
 242 //      W1 = ----       W2 = ---        W3 = -------
 243 //             2              x               x + 1
 244 //            x
 245 //
 246 // Then the coefficient matrix C is
 247 //
 248 //              coeff(W1,x,0)   coeff(W2,x,0)   coeff(W3,x,0)
 249 //
 250 //       C =    coeff(W1,x,1)   coeff(W2,x,1)   coeff(W3,x,1)
 251 //
 252 //              coeff(W1,x,2)   coeff(W2,x,2)   coeff(W3,x,2)
 253 //
 254 // It follows that
 255 //
 256 //       coeff(B,x,0)     Y1
 257 //
 258 //       coeff(B,x,1) = C Y2
 259 //
 260 //       coeff(B,x,2) =   Y3
 261 //
 262 // Hence
 263 //
 264 //       Y1       coeff(B,x,0)
 265 //             -1
 266 //       Y2 = C   coeff(B,x,1)
 267 //
 268 //       Y3       coeff(B,x,2)
 269
 270 void
 271 expand_get_C(void)
 272 {
 273         int h, i, j, n;
 274         U **a;
 275         h = tos;
 276         if (car(A) == symbol(MULTIPLY)) {
 277                 p1 = cdr(A);
 278                 while (iscons(p1)) {
 279                         F = car(p1);
 280                         expand_get_CF();
 281                         p1 = cdr(p1);
 282                 }
 283         } else {
 284                 F = A;
 285                 expand_get_CF();
 286         }
 287         n = tos - h;
 288         if (n == 1) {
 289                 C = pop();
 290                 return;
 291         }
 292         C = alloc_tensor(n * n);
 293         C->u.tensor->ndim = 2;
 294         C->u.tensor->dim[0] = n;
 295         C->u.tensor->dim[1] = n;
 296         a = stack + h;
 297         for (i = 0; i < n; i++) {
 298                 for (j = 0; j < n; j++) {
 299                         push(a[j]);
 300                         push(X);
 301                         push_integer(i);
 302                         power();
 303                         divide();
 304                         push(X);
 305                         filter();
 306                         C->u.tensor->elem[n * i + j] = pop();
 307                 }
 308         }
 309         tos -= n;
 310 }
 311
 312 // The following table shows the push order for simple roots, repeated roots,
 313 // and inrreducible factors.
 314 //
 315 //  Factor F        Push 1st        Push 2nd         Push 3rd      Push 4th
 316 //
 317 //
 318 //                   A
 319 //  x               ---
 320 //                   x
 321 //
 322 //
 323 //   2               A               A
 324 //  x               ----            ---
 325 //                    2              x
 326 //                   x
 327 //
 328 //
 329 //                     A
 330 //  x + 1           -------
 331 //                   x + 1
 332 //
 333 //
 334 //         2            A              A
 335 //  (x + 1)         ----------      -------
 336 //                          2        x + 1
 337 //                   (x + 1)
 338 //
 339 //
 340 //   2                   A               Ax
 341 //  x  + x + 1      ------------    ------------
 342 //                    2               2
 343 //                   x  + x + 1      x  + x + 1
 344 //
 345 //
 346 //    2         2          A              Ax              A             Ax
 347 //  (x  + x + 1)    --------------- ---------------  ------------  ------------
 348 //                     2         2     2         2     2             2
 349 //                   (x  + x + 1)    (x  + x + 1)     x  + x + 1    x  + x + 1
 350 //
 351 //
 352 // For T = A/F and F = P^N we have
 353 //
 354 //
 355 //      Factor F          Push 1st    Push 2nd    Push 3rd    Push 4th
 356 //
 357 //      x                 T
 358 //
 359 //       2
 360 //      x                 T           TP
 361 //
 362 //
 363 //      x + 1             T
 364 //
 365 //             2
 366 //      (x + 1)           T           TP
 367 //
 368 //       2
 369 //      x  + x + 1        T           TX
 370 //
 371 //        2         2
 372 //      (x  + x + 1)      T           TX          TP          TPX
 373 //
 374 //
 375 // Hence we want to push in the order
 376 //
 377 //      T * (P ^ i) * (X ^ j)
 378 //
 379 // for all i, j such that
 380 //
 381 //      i = 0, 1, ..., N - 1
 382 //
 383 //      j = 0, 1, ..., deg(P) - 1
 384 //
 385 // where index j runs first.
 386
 387 void
 388 expand_get_CF(void)
 389 {       int d, i, j, n;
 390         if (!find(F, X))
 391                 return;
 392         trivial_divide();
 393         if (car(F) == symbol(POWER)) {
 394                 push(caddr(F));
 395                 n = pop_integer();
 396                 P = cadr(F);
 397         } else {
 398                 n = 1;
 399                 P = F;
 400         }
 401         push(P);
 402         push(X);
 403         degree();
 404         d = pop_integer();
 405         for (i = 0; i < n; i++) {
 406                 for (j = 0; j < d; j++) {
 407                         push(T);
 408                         push(P);
 409                         push_integer(i);
 410                         power();
 411                         multiply();
 412                         push(X);
 413                         push_integer(j);
 414                         power();
 415                         multiply();
 416                 }
 417         }
 418 }
 419
 420 // Returns T = A/F where F is a factor of A.
 421
 422 void
 423 trivial_divide(void)
 424 {
 425         int h;
 426         if (car(A) == symbol(MULTIPLY)) {
 427                 h = tos;
 428                 p0 = cdr(A);
 429                 while (iscons(p0)) {
 430                         if (!equal(car(p0), F)) {
 431                                 push(car(p0));
 432                                 eval(); // force expansion of (x+1)^2, f.e.
 433                         }
 434                         p0 = cdr(p0);
 435                 }
 436                 multiply_all(tos - h);
 437         } else
 438                 push_integer(1);
 439         T = pop();
 440 }
 441
 442 // Returns the expansion coefficient vector B.
 443
 444 void
 445 expand_get_B(void)
 446 {
 447         int i, n;
 448         if (!istensor(C))
 449                 return;
 450         n = C->u.tensor->dim[0];
 451         T = alloc_tensor(n);
 452         T->u.tensor->ndim = 1;
 453         T->u.tensor->dim[0] = n;
 454         for (i = 0; i < n; i++) {
 455                 push(B);
 456                 push(X);
 457                 push_integer(i);
 458                 power();
 459                 divide();
 460                 push(X);
 461                 filter();
 462                 T->u.tensor->elem[i] = pop();
 463         }
 464         B = T;
 465 }
 466
 467 // Returns the expansion fractions in A.
 468
 469 void
 470 expand_get_A(void)
 471 {
 472         int h, i, n;
 473         if (!istensor(C)) {
 474                 push(A);
 475                 reciprocate();
 476                 A = pop();
 477                 return;
 478         }
 479         h = tos;
 480         if (car(A) == symbol(MULTIPLY)) {
 481                 T = cdr(A);
 482                 while (iscons(T)) {
 483                         F = car(T);
 484                         expand_get_AF();
 485                         T = cdr(T);
 486                 }
 487         } else {
 488                 F = A;
 489                 expand_get_AF();
 490         }
 491         n = tos - h;
 492         T = alloc_tensor(n);
 493         T->u.tensor->ndim = 1;
 494         T->u.tensor->dim[0] = n;
 495         for (i = 0; i < n; i++)
 496                 T->u.tensor->elem[i] = stack[h + i];
 497         tos = h;
 498         A = T;
 499 }
 500
 501 void
 502 expand_get_AF(void)
 503 {       int d, i, j, n = 1;
 504         if (!find(F, X))
 505                 return;
 506         if (car(F) == symbol(POWER)) {
 507                 push(caddr(F));
 508                 n = pop_integer();
 509                 F = cadr(F);
 510         }
 511         push(F);
 512         push(X);
 513         degree();
 514         d = pop_integer();
 515         for (i = n; i > 0; i--) {
 516                 for (j = 0; j < d; j++) {
 517                         push(F);
 518                         push_integer(i);
 519                         power();
 520                         reciprocate();
 521                         push(X);
 522                         push_integer(j);
 523                         power();
 524                         multiply();
 525                 }
 526         }
 527 }
 528
 529 #if SELFTEST
 530
 531 static char *s[] = {
 532
 533         // general cases
 534
 535         "expand(1/(x+1)/(x+2))",
 536         "1/(x+1)-1/(x+2)",
 537
 538         "expand((2x^3-x+2)/(x^2-2x+1))",
 539         "4+2*x+5/(x-1)+3/(x^2-2*x+1)",
 540
 541         "expand(1/x^2/(x-1))",
 542         "-1/(x^2)-1/x+1/(x-1)",
 543
 544         "p=5s+2",
 545         "",
 546
 547         "q=(s+1)(s+2)^2",
 548         "",
 549
 550         "expand(p/q)",
 551         "-3/(s+1)+3/(s+2)+8/(s^2+4*s+4)",
 552
 553         // ensure denominators are expanded (result seems preferable that way)
 554
 555         "q=(x-1)(x-2)^3",
 556         "",
 557
 558         "expand(1/q)",
 559         "1/(x^3-6*x^2+12*x-8)+1/(x-2)-1/(x-1)-1/(x^2-4*x+4)",
 560
 561         // fractional poles
 562
 563         "expand(1/(x+1/2)/(x+1/3))",
 564         "-12/(2*x+1)+18/(3*x+1)",
 565
 566         // expand tensor
 567
 568         "f=1/(x+1)/(x+2)",
 569         "",
 570
 571         "g=1/(x+1)-1/(x+2)",
 572         "",
 573
 574         "expand(((f,f),(f,f)))-((g,g),(g,g))",
 575         "((0,0),(0,0))",
 576
 577         // denominator normalized?
 578
 579         "expand(1/(1+1/x))",
 580         "1-1/(x+1)",
 581
 582         // poles at zero
 583
 584         "expand(1/x/(x+1))",
 585         "1/x-1/(x+1)",
 586
 587         "expand(1/x^2/(x+1))",
 588         "x^(-2)-1/x+1/(x+1)",
 589
 590         // other corner cases
 591
 592         "expand(1/x)",
 593         "1/x",
 594
 595         "expand(1/x^2)",
 596         "x^(-2)",
 597
 598         "expand(1/(x^2-4x+4))",
 599         "1/(x^2-4*x+4)",
 600 };
 601
 602 void
 603 test_expand(void)
 604 {
 605         test(__FILE__, s, sizeof s / sizeof (char *));
 606 }
 607
 608 #endif