multiply.cpp

   1 // Symbolic multiplication
   2
   3 #include "stdafx.h"
   4 #include "defs.h"
   5
   6 extern void append(void);
   7 static void parse_p1(void);
   8 static void parse_p2(void);
   9 static void __normalize_radical_factors(int);
  10
  11 void
  12 multiply(void)
  13 {
  14         if (esc_flag)
  15                 stop("escape key stop");
  16         if (isnum(stack[tos - 2]) && isnum(stack[tos - 1]))
  17                 multiply_numbers();
  18         else {
  19                 save();
  20                 yymultiply();
  21                 restore();
  22         }
  23 }
  24
  25 void
  26 yymultiply(void)
  27 {
  28         int h, i, n;
  29
  30         // pop operands
  31
  32         p2 = pop();
  33         p1 = pop();
  34
  35         h = tos;
  36
  37         // is either operand zero?
  38
  39         if (iszero(p1) || iszero(p2)) {
  40                 push(zero);
  41                 return;
  42         }
  43
  44         // is either operand a sum?
  45
  46         if (expanding && isadd(p1)) {
  47                 p1 = cdr(p1);
  48                 push(zero);
  49                 while (iscons(p1)) {
  50                         push(car(p1));
  51                         push(p2);
  52                         multiply();
  53                         add();
  54                         p1 = cdr(p1);
  55                 }
  56                 return;
  57         }
  58
  59         if (expanding && isadd(p2)) {
  60                 p2 = cdr(p2);
  61                 push(zero);
  62                 while (iscons(p2)) {
  63                         push(p1);
  64                         push(car(p2));
  65                         multiply();
  66                         add();
  67                         p2 = cdr(p2);
  68                 }
  69                 return;
  70         }
  71
  72         // scalar times tensor?
  73
  74         if (!istensor(p1) && istensor(p2)) {
  75                 push(p1);
  76                 push(p2);
  77                 scalar_times_tensor();
  78                 return;
  79         }
  80
  81         // tensor times scalar?
  82
  83         if (istensor(p1) && !istensor(p2)) {
  84                 push(p1);
  85                 push(p2);
  86                 tensor_times_scalar();
  87                 return;
  88         }
  89
  90         // adjust operands
  91
  92         if (car(p1) == symbol(MULTIPLY))
  93                 p1 = cdr(p1);
  94         else {
  95                 push(p1);
  96                 list(1);
  97                 p1 = pop();
  98         }
  99
 100         if (car(p2) == symbol(MULTIPLY))
 101                 p2 = cdr(p2);
 102         else {
 103                 push(p2);
 104                 list(1);
 105                 p2 = pop();
 106         }
 107
 108         // handle numerical coefficients
 109
 110         if (isnum(car(p1)) && isnum(car(p2))) {
 111                 push(car(p1));
 112                 push(car(p2));
 113                 multiply_numbers();
 114                 p1 = cdr(p1);
 115                 p2 = cdr(p2);
 116         } else if (isnum(car(p1))) {
 117                 push(car(p1));
 118                 p1 = cdr(p1);
 119         } else if (isnum(car(p2))) {
 120                 push(car(p2));
 121                 p2 = cdr(p2);
 122         } else
 123                 push(one);
 124
 125         parse_p1();
 126         parse_p2();
 127
 128         while (iscons(p1) && iscons(p2)) {
 129
 130 //              if (car(p1)->gamma && car(p2)->gamma) {
 131 //                      combine_gammas(h);
 132 //                      p1 = cdr(p1);
 133 //                      p2 = cdr(p2);
 134 //                      parse_p1();
 135 //                      parse_p2();
 136 //                      continue;
 137 //              }
 138
 139                 if (caar(p1) == symbol(OPERATOR) && caar(p2) == symbol(OPERATOR)) {
 140                         push_symbol(OPERATOR);
 141                         push(cdar(p1));
 142                         push(cdar(p2));
 143                         append();
 144                         cons();
 145                         p1 = cdr(p1);
 146                         p2 = cdr(p2);
 147                         parse_p1();
 148                         parse_p2();
 149                         continue;
 150                 }
 151
 152                 switch (cmp_expr(p3, p4)) {
 153                 case -1:
 154                         push(car(p1));
 155                         p1 = cdr(p1);
 156                         parse_p1();
 157                         break;
 158                 case 1:
 159                         push(car(p2));
 160                         p2 = cdr(p2);
 161                         parse_p2();
 162                         break;
 163                 case 0:
 164                         combine_factors(h);
 165                         p1 = cdr(p1);
 166                         p2 = cdr(p2);
 167                         parse_p1();
 168                         parse_p2();
 169                         break;
 170                 default:
 171                         stop("internal error 2");
 172                         break;
 173                 }
 174         }
 175
 176         // push remaining factors, if any
 177
 178         while (iscons(p1)) {
 179                 push(car(p1));
 180                 p1 = cdr(p1);
 181         }
 182
 183         while (iscons(p2)) {
 184                 push(car(p2));
 185                 p2 = cdr(p2);
 186         }
 187
 188         // normalize radical factors
 189
 190         // example: 2*2(-1/2) -> 2^(1/2)
 191
 192         // must be done after merge because merge may produce radical
 193
 194         // example: 2^(1/2-a)*2^a -> 2^(1/2)
 195
 196         __normalize_radical_factors(h);
 197
 198         // this hack should not be necessary, unless power returns a multiply
 199
 200         //for (i = h; i < tos; i++) {
 201         //      if (car(stack[i]) == symbol(MULTIPLY)) {
 202         //              multiply_all(tos - h);
 203         //              return;
 204         //      }
 205         //}
 206
 207         if (expanding) {
 208                 for (i = h; i < tos; i++) {
 209                         if (isadd(stack[i])) {
 210                                 multiply_all(tos - h);
 211                                 return;
 212                         }
 213                 }
 214         }
 215
 216         // n is the number of result factors on the stack
 217
 218         n = tos - h;
 219
 220         if (n == 1)
 221                 return;
 222
 223         // discard integer 1
 224
 225         if (isrational(stack[h]) && equaln(stack[h], 1)) {
 226                 if (n == 2) {
 227                         p7 = pop();
 228                         pop();
 229                         push(p7);
 230                 } else {
 231                         stack[h] = symbol(MULTIPLY);
 232                         list(n);
 233                 }
 234                 return;
 235         }
 236
 237         list(n);
 238         p7 = pop();
 239         push_symbol(MULTIPLY);
 240         push(p7);
 241         cons();
 242 }
 243
 244 // Decompose a factor into base and power.
 245 //
 246 // input:       car(p1)         factor
 247 //
 248 // output:      p3              factor's base
 249 //
 250 //              p5              factor's power (possibly 1)
 251
 252 static void
 253 parse_p1(void)
 254 {
 255         p3 = car(p1);
 256         p5 = one;
 257         if (car(p3) == symbol(POWER)) {
 258                 p5 = caddr(p3);
 259                 p3 = cadr(p3);
 260         }
 261 }
 262
 263 // Decompose a factor into base and power.
 264 //
 265 // input:       car(p2)         factor
 266 //
 267 // output:      p4              factor's base
 268 //
 269 //              p6              factor's power (possibly 1)
 270
 271 static void
 272 parse_p2(void)
 273 {
 274         p4 = car(p2);
 275         p6 = one;
 276         if (car(p4) == symbol(POWER)) {
 277                 p6 = caddr(p4);
 278                 p4 = cadr(p4);
 279         }
 280 }
 281
 282 void
 283 combine_factors(int h)
 284 {
 285         push(p4);
 286         push(p5);
 287         push(p6);
 288         add();
 289         power();
 290         p7 = pop();
 291         if (isnum(p7)) {
 292                 push(stack[h]);
 293                 push(p7);
 294                 multiply_numbers();
 295                 stack[h] = pop();
 296         } else if (car(p7) == symbol(MULTIPLY)) {
 297                 // power can return number * factor (i.e. -1 * i)
 298                 if (isnum(cadr(p7)) && cdddr(p7) == symbol(NIL)) {
 299                         push(stack[h]);
 300                         push(cadr(p7));
 301                         multiply_numbers();
 302                         stack[h] = pop();
 303                         push(caddr(p7));
 304                 } else
 305                         push(p7);
 306         } else
 307                 push(p7);
 308 }
 309
 310 int gp[17][17] = {
 311         {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
 312         {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
 313         {0,0,1,-6,-7,-8,-3,-4,-5,13,14,15,-16,9,10,11,-12},
 314         {0,0,6,-1,-11,10,-2,-15,14,12,-5,4,-9,16,-8,7,-13},
 315         {0,0,7,11,-1,-9,15,-2,-13,5,12,-3,-10,8,16,-6,-14},
 316         {0,0,8,-10,9,-1,-14,13,-2,-4,3,12,-11,-7,6,16,-15},
 317         {0,0,3,2,15,-14,1,11,-10,16,-8,7,13,12,-5,4,9},
 318         {0,0,4,-15,2,13,-11,1,9,8,16,-6,14,5,12,-3,10},
 319         {0,0,5,14,-13,2,10,-9,1,-7,6,16,15,-4,3,12,11},
 320         {0,0,13,12,-5,4,16,-8,7,-1,-11,10,-3,-2,-15,14,-6},
 321         {0,0,14,5,12,-3,8,16,-6,11,-1,-9,-4,15,-2,-13,-7},
 322         {0,0,15,-4,3,12,-7,6,16,-10,9,-1,-5,-14,13,-2,-8},
 323         {0,0,16,-9,-10,-11,-13,-14,-15,-3,-4,-5,1,-6,-7,-8,2},
 324         {0,0,9,-16,8,-7,-12,5,-4,-2,-15,14,6,-1,-11,10,3},
 325         {0,0,10,-8,-16,6,-5,-12,3,15,-2,-13,7,11,-1,-9,4},
 326         {0,0,11,7,-6,-16,4,-3,-12,-14,13,-2,8,-10,9,-1,5},
 327         {0,0,12,13,14,15,9,10,11,-6,-7,-8,-2,-3,-4,-5,-1}
 328 };
 329
 330 #if 0
 331
 332 static void
 333 combine_gammas(int h)
 334 {
 335         int n;
 336         n = gp[(int) p1->gamma][(int) p2->gamma];
 337         if (n < 0) {
 338                 n = -n;
 339                 push(stack[h]);
 340                 negate();
 341                 stack[h] = pop();
 342         }
 343         if (n > 1)
 344                 push(_gamma[n]);
 345 }
 346
 347 #endif
 348
 349 void
 350 multiply_noexpand(void)
 351 {
 352         int x;
 353         x = expanding;
 354         expanding = 0;
 355         multiply();
 356         expanding = x;
 357 }
 358
 359 // multiply n factors on stack
 360
 361 void
 362 multiply_all(int n)
 363 {
 364         int h, i;
 365         if (n == 1)
 366                 return;
 367         if (n == 0) {
 368                 push(one);
 369                 return;
 370         }
 371         h = tos - n;
 372         push(stack[h]);
 373         for (i = 1; i < n; i++) {
 374                 push(stack[h + i]);
 375                 multiply();
 376         }
 377         stack[h] = pop();
 378         tos = h + 1;
 379 }
 380
 381 void
 382 multiply_all_noexpand(int n)
 383 {
 384         int x;
 385         x = expanding;
 386         expanding = 0;
 387         multiply_all(n);
 388         expanding = x;
 389 }
 390
 391 //-----------------------------------------------------------------------------
 392 //
 393 //      Symbolic division
 394 //
 395 //      Input:          Dividend and divisor on stack
 396 //
 397 //      Output:         Quotient on stack
 398 //
 399 //-----------------------------------------------------------------------------
 400
 401 void
 402 divide(void)
 403 {
 404         if (isnum(stack[tos - 2]) && isnum(stack[tos - 1]))
 405                 divide_numbers();
 406         else {
 407                 inverse();
 408                 multiply();
 409         }
 410 }
 411
 412 void
 413 inverse(void)
 414 {
 415         if (isnum(stack[tos - 1]))
 416                 invert_number();
 417         else {
 418                 push_integer(-1);
 419                 power();
 420         }
 421 }
 422
 423 void
 424 reciprocate(void)
 425 {
 426         if (isnum(stack[tos - 1]))
 427                 invert_number();
 428         else {
 429                 push_integer(-1);
 430                 power();
 431         }
 432 }
 433
 434 void
 435 negate(void)
 436 {
 437         if (isnum(stack[tos - 1]))
 438                 negate_number();
 439         else {
 440                 push_integer(-1);
 441                 multiply();
 442         }
 443 }
 444
 445 void
 446 negate_expand(void)
 447 {
 448         int x;
 449         x = expanding;
 450         expanding = 1;
 451         negate();
 452         expanding = x;
 453 }
 454
 455 void
 456 negate_noexpand(void)
 457 {
 458         int x;
 459         x = expanding;
 460         expanding = 0;
 461         negate();
 462         expanding = x;
 463 }
 464
 465 //-----------------------------------------------------------------------------
 466 //
 467 //      Normalize radical factors
 468 //
 469 //      Input:          stack[h]        Coefficient factor, possibly 1
 470 //
 471 //                      stack[h + 1]    Second factor
 472 //
 473 //                      stack[tos - 1]  Last factor
 474 //
 475 //      Output:         Reduced coefficent and normalized radicals (maybe)
 476 //
 477 //      Example:        2*2^(-1/2) -> 2^(1/2)
 478 //
 479 //      (power number number) is guaranteed to have the following properties:
 480 //
 481 //      1. Base is an integer
 482 //
 483 //      2. Absolute value of exponent < 1
 484 //
 485 //      These properties are assured by the power function.
 486 //
 487 //-----------------------------------------------------------------------------
 488
 489 #define A p1
 490 #define B p2
 491
 492 #define BASE p3
 493 #define EXPO p4
 494
 495 #define TMP p5
 496
 497 static int __is_radical_number(U *);
 498
 499 static void
 500 __normalize_radical_factors(int h)
 501 {
 502         int i;
 503
 504         // if coeff is 1 or floating then don't bother
 505
 506         if (isplusone(stack[h]) || isminusone(stack[h]) || isdouble(stack[h]))
 507                 return;
 508
 509         // if no radicals then don't bother
 510
 511         for (i = h + 1; i < tos; i++)
 512                 if (__is_radical_number(stack[i]))
 513                         break;
 514
 515         if (i == tos)
 516                 return;
 517
 518         // ok, try to simplify
 519
 520         save();
 521
 522         // numerator
 523
 524         push(stack[h]);
 525         mp_numerator();
 526         A = pop();
 527
 528         for (i = h + 1; i < tos; i++) {
 529
 530                 if (isplusone(A) || isminusone(A))
 531                         break;
 532
 533                 if (!__is_radical_number(stack[i]))
 534                         continue;
 535
 536                 BASE = cadr(stack[i]);
 537                 EXPO = caddr(stack[i]);
 538
 539                 // exponent must be negative
 540
 541                 if (!isnegativenumber(EXPO))
 542                         continue;
 543
 544                 // numerator divisible by BASE?
 545
 546                 push(A);
 547                 push(BASE);
 548                 divide();
 549
 550                 TMP = pop();
 551
 552                 if (!isinteger(TMP))
 553                         continue;
 554
 555                 // reduce numerator
 556
 557                 A = TMP;
 558
 559                 // invert radical
 560
 561                 push_symbol(POWER);
 562                 push(BASE);
 563                 push(one);
 564                 push(EXPO);
 565                 add();
 566                 list(3);
 567                 stack[i] = pop();
 568         }
 569
 570         // denominator
 571
 572         push(stack[h]);
 573         mp_denominator();
 574         B = pop();
 575
 576         for (i = h + 1; i < tos; i++) {
 577
 578                 if (isplusone(B))
 579                         break;
 580
 581                 if (!__is_radical_number(stack[i]))
 582                         continue;
 583
 584                 BASE = cadr(stack[i]);
 585                 EXPO = caddr(stack[i]);
 586
 587                 // exponent must be positive
 588
 589                 if (isnegativenumber(EXPO))
 590                         continue;
 591
 592                 // denominator divisible by BASE?
 593
 594                 push(B);
 595                 push(BASE);
 596                 divide();
 597
 598                 TMP = pop();
 599
 600                 if (!isinteger(TMP))
 601                         continue;
 602
 603                 // reduce denominator
 604
 605                 B = TMP;
 606
 607                 // invert radical
 608
 609                 push_symbol(POWER);
 610                 push(BASE);
 611                 push(EXPO);
 612                 push(one);
 613                 subtract();
 614                 list(3);
 615                 stack[i] = pop();
 616         }
 617
 618         // reconstitute the coefficient
 619
 620         push(A);
 621         push(B);
 622         divide();
 623         stack[h] = pop();
 624
 625         restore();
 626 }
 627
 628 // don't include i
 629
 630 static int
 631 __is_radical_number(U *p)
 632 {
 633         // don't use i
 634
 635         if (car(p) == symbol(POWER) && isnum(cadr(p)) && isnum(caddr(p)) && !isminusone(cadr(p)))
 636                 return 1;
 637         else
 638                 return 0;
 639 }
 640
 641 //-----------------------------------------------------------------------------
 642 //
 643 //      > a*hilbert(2)
 644 //      ((a,1/2*a),(1/2*a,1/3*a))
 645 //
 646 //      Note that "a" is presumed to be a scalar. Is this correct?
 647 //
 648 //      Yes, because "*" has no meaning if "a" is a tensor.
 649 //      To multiply tensors, "dot" or "outer" should be used.
 650 //
 651 //      > dot(a,hilbert(2))
 652 //      dot(a,((1,1/2),(1/2,1/3)))
 653 //
 654 //      In this case "a" could be a scalar or tensor so the result is not
 655 //      expanded.
 656 //
 657 //-----------------------------------------------------------------------------
 658
 659 #if SELFTEST
 660
 661 static char *s[] = {
 662
 663         "0*a",
 664         "0",
 665
 666         "a*0",
 667         "0",
 668
 669         "1*a",
 670         "a",
 671
 672         "a*1",
 673         "a",
 674
 675         "a*a",
 676         "a^2",
 677
 678         "a^2*a",
 679         "a^3",
 680
 681         "a*a^2",
 682         "a^3",
 683
 684         "a^2*a^2",
 685         "a^4",
 686
 687         "2^a*2^(3-a)",          // symbolic exponents cancel
 688         "8",
 689
 690         "sqrt(2)/2",
 691         "2^(-1/2)",
 692
 693         "2/sqrt(2)",
 694         "2^(1/2)",
 695
 696         "-sqrt(2)/2",
 697         "-1/(2^(1/2))",
 698
 699         "2^(1/2-a)*2^a/10",
 700         "1/(5*2^(1/2))",
 701
 702         "i/4",
 703         "1/4*i",
 704
 705         "1/(4 i)",
 706         "-1/4*i",
 707
 708         // ensure 1.0 is not discarded
 709
 710         "1.0 pi 1/2",
 711         "0.5*pi",
 712
 713         "1.0 1/2 pi",
 714         "0.5*pi",
 715 };
 716
 717 void
 718 test_multiply(void)
 719 {
 720         test(__FILE__, s, sizeof s / sizeof (char *));
 721 }
 722
 723 #endif