power.cpp

   1 /* Power function
   2
   3         Input:          push    Base
   4
   5                         push    Exponent
   6
   7         Output:         Result on stack
   8 */
   9
  10 #include "stdafx.h"
  11 #include "defs.h"
  12
  13 void
  14 eval_power(void)
  15 {
  16         push(cadr(p1));
  17         eval();
  18         push(caddr(p1));
  19         eval();
  20         power();
  21 }
  22
  23 void
  24 power(void)
  25 {
  26         save();
  27         yypower();
  28         restore();
  29 }
  30
  31 void
  32 yypower(void)
  33 {
  34         int n;
  35
  36         p2 = pop();
  37         p1 = pop();
  38
  39         // both base and exponent are rational numbers?
  40
  41         if (isrational(p1) && isrational(p2)) {
  42                 push(p1);
  43                 push(p2);
  44                 qpow();
  45                 return;
  46         }
  47
  48         // both base and exponent are either rational or double?
  49
  50         if (isnum(p1) && isnum(p2)) {
  51                 push(p1);
  52                 push(p2);
  53                 dpow();
  54                 return;
  55         }
  56
  57         if (istensor(p1)) {
  58                 power_tensor();
  59                 return;
  60         }
  61
  62         if (p1 == symbol(E) && car(p2) == symbol(LOG)) {
  63                 push(cadr(p2));
  64                 return;
  65         }
  66
  67         if (p1 == symbol(E) && isdouble(p2)) {
  68                 push_double(exp(p2->u.d));
  69                 return;
  70         }
  71
  72         //      1 ^ a           ->      1
  73
  74         //      a ^ 0           ->      1
  75
  76         if (equal(p1, one) || iszero(p2)) {
  77                 push(one);
  78                 return;
  79         }
  80
  81         //      a ^ 1           ->      a
  82
  83         if (equal(p2, one)) {
  84                 push(p1);
  85                 return;
  86         }
  87
  88         //      (a * b) ^ c     ->      (a ^ c) * (b ^ c)
  89
  90         if (car(p1) == symbol(MULTIPLY)) {
  91                 p1 = cdr(p1);
  92                 push(car(p1));
  93                 push(p2);
  94                 power();
  95                 p1 = cdr(p1);
  96                 while (iscons(p1)) {
  97                         push(car(p1));
  98                         push(p2);
  99                         power();
 100                         multiply();
 101                         p1 = cdr(p1);
 102                 }
 103                 return;
 104         }
 105
 106         //      (a ^ b) ^ c     ->      a ^ (b * c)
 107
 108         if (car(p1) == symbol(POWER)) {
 109                 push(cadr(p1));
 110                 push(caddr(p1));
 111                 push(p2);
 112                 multiply();
 113                 power();
 114                 return;
 115         }
 116
 117         //      (a + b) ^ n     ->      (a + b) * (a + b) ...
 118
 119         if (expanding && isadd(p1) && isnum(p2)) {
 120                 push(p2);
 121                 n = pop_integer();
 122                 if (n > 1) {
 123                         power_sum(n);
 124                         return;
 125                 }
 126         }
 127
 128         //      sin(x) ^ 2n -> (1 - cos(x) ^ 2) ^ n
 129
 130         if (trigmode == 1 && car(p1) == symbol(SIN) && iseveninteger(p2)) {
 131                 push_integer(1);
 132                 push(cadr(p1));
 133                 cosine();
 134                 push_integer(2);
 135                 power();
 136                 subtract();
 137                 push(p2);
 138                 push_rational(1, 2);
 139                 multiply();
 140                 power();
 141                 return;
 142         }
 143
 144         //      cos(x) ^ 2n -> (1 - sin(x) ^ 2) ^ n
 145
 146         if (trigmode == 2 && car(p1) == symbol(COS) && iseveninteger(p2)) {
 147                 push_integer(1);
 148                 push(cadr(p1));
 149                 sine();
 150                 push_integer(2);
 151                 power();
 152                 subtract();
 153                 push(p2);
 154                 push_rational(1, 2);
 155                 multiply();
 156                 power();
 157                 return;
 158         }
 159
 160         // complex number? (just number, not expression)
 161
 162         if (iscomplexnumber(p1)) {
 163
 164                 // integer power?
 165
 166                 // n will be negative here, positive n already handled
 167
 168                 if (isinteger(p2)) {
 169
 170                         //               /        \  n
 171                         //         -n   |  a - ib  |
 172                         // (a + ib)   = | -------- |
 173                         //              |   2   2  |
 174                         //               \ a + b  /
 175
 176                         push(p1);
 177                         conjugate();
 178                         p3 = pop();
 179                         push(p3);
 180                         push(p3);
 181                         push(p1);
 182                         multiply();
 183                         divide();
 184                         push(p2);
 185                         negate();
 186                         power();
 187                         return;
 188                 }
 189
 190                 // noninteger or floating power?
 191
 192                 if (isnum(p2)) {
 193
 194 #if 1                   // use polar form
 195                         push(p1);
 196                         mag();
 197                         push(p2);
 198                         power();
 199                         push_integer(-1);
 200                         push(p1);
 201                         arg();
 202                         push(p2);
 203                         multiply();
 204                         push(symbol(PI));
 205                         divide();
 206                         power();
 207                         multiply();
 208
 209 #else                   // use exponential form
 210                         push(p1);
 211                         mag();
 212                         push(p2);
 213                         power();
 214                         push(symbol(E));
 215                         push(p1);
 216                         arg();
 217                         push(p2);
 218                         multiply();
 219                         push(imaginaryunit);
 220                         multiply();
 221                         power();
 222                         multiply();
 223 #endif
 224                         return;
 225                 }
 226         }
 227
 228         if (simplify_polar())
 229                 return;
 230
 231         push_symbol(POWER);
 232         push(p1);
 233         push(p2);
 234         list(3);
 235 }
 236
 237 //-----------------------------------------------------------------------------
 238 //
 239 //      Compute the power of a sum
 240 //
 241 //      Input:          p1      sum
 242 //
 243 //                      n       exponent
 244 //
 245 //      Output:         Result on stack
 246 //
 247 //      Note:
 248 //
 249 //      Uses the multinomial series (see Math World)
 250 //
 251 //                          n              n!          n1   n2       nk
 252 //      (a1 + a2 + ... + ak)  = sum (--------------- a1   a2   ... ak  )
 253 //                                   n1! n2! ... nk!
 254 //
 255 //      The sum is over all n1 ... nk such that n1 + n2 + ... + nk = n.
 256 //
 257 //-----------------------------------------------------------------------------
 258
 259 // first index is the term number 0..k-1, second index is the exponent 0..n
 260
 261 #define A(i, j) frame[(i) * (n + 1) + (j)]
 262
 263 void
 264 power_sum(int n)
 265 {
 266         int *a, i, j, k;
 267
 268         // number of terms in the sum
 269
 270         k = length(p1) - 1;
 271
 272         // local frame
 273
 274         push_frame(k * (n + 1));
 275
 276         // array of powers
 277
 278         p1 = cdr(p1);
 279         for (i = 0; i < k; i++) {
 280                 for (j = 0; j <= n; j++) {
 281                         push(car(p1));
 282                         push_integer(j);
 283                         power();
 284                         A(i, j) = pop();
 285                 }
 286                 p1 = cdr(p1);
 287         }
 288
 289         push_integer(n);
 290         factorial();
 291         p1 = pop();
 292
 293         a = (int *) malloc(k * sizeof (int));
 294
 295         if (a == NULL)
 296                 stop("malloc failure");
 297
 298         push(zero);
 299
 300         multinomial_sum(k, n, a, 0, n);
 301
 302         free(a);
 303
 304         pop_frame(k * (n + 1));
 305 }
 306
 307 //-----------------------------------------------------------------------------
 308 //
 309 //      Compute multinomial sum
 310 //
 311 //      Input:          k       number of factors
 312 //
 313 //                      n       overall exponent
 314 //
 315 //                      a       partition array
 316 //
 317 //                      i       partition array index
 318 //
 319 //                      m       partition remainder
 320 //
 321 //                      p1      n!
 322 //
 323 //                      A       factor array
 324 //
 325 //      Output:         Result on stack
 326 //
 327 //      Note:
 328 //
 329 //      Uses recursive descent to fill the partition array.
 330 //
 331 //-----------------------------------------------------------------------------
 332
 333 void
 334 multinomial_sum(int k, int n, int *a, int i, int m)
 335 {
 336         int j;
 337
 338         if (i < k - 1) {
 339                 for (j = 0; j <= m; j++) {
 340                         a[i] = j;
 341                         multinomial_sum(k, n, a, i + 1, m - j);
 342                 }
 343                 return;
 344         }
 345
 346         a[i] = m;
 347
 348         // coefficient
 349
 350         push(p1);
 351
 352         for (j = 0; j < k; j++) {
 353                 push_integer(a[j]);
 354                 factorial();
 355                 divide();
 356         }
 357
 358         // factors
 359
 360         for (j = 0; j < k; j++) {
 361                 push(A(j, a[j]));
 362                 multiply();
 363         }
 364
 365         add();
 366 }
 367
 368 // exp(n/2 i pi) ?
 369
 370 // p2 is the exponent expression
 371
 372 // clobbers p3
 373
 374 int
 375 simplify_polar(void)
 376 {
 377         int n;
 378
 379         n = isquarterturn(p2);
 380
 381         switch(n) {
 382         case 0:
 383                 break;
 384         case 1:
 385                 push_integer(1);
 386                 return 1;
 387         case 2:
 388                 push_integer(-1);
 389                 return 1;
 390         case 3:
 391                 push(imaginaryunit);
 392                 return 1;
 393         case 4:
 394                 push(imaginaryunit);
 395                 negate();
 396                 return 1;
 397         }
 398
 399         if (car(p2) == symbol(ADD)) {
 400                 p3 = cdr(p2);
 401                 while (iscons(p3)) {
 402                         n = isquarterturn(car(p3));
 403                         if (n)
 404                                 break;
 405                         p3 = cdr(p3);
 406                 }
 407                 switch (n) {
 408                 case 0:
 409                         return 0;
 410                 case 1:
 411                         push_integer(1);
 412                         break;
 413                 case 2:
 414                         push_integer(-1);
 415                         break;
 416                 case 3:
 417                         push(imaginaryunit);
 418                         break;
 419                 case 4:
 420                         push(imaginaryunit);
 421                         negate();
 422                         break;
 423                 }
 424                 push(p2);
 425                 push(car(p3));
 426                 subtract();
 427                 exponential();
 428                 multiply();
 429                 return 1;
 430         }
 431
 432         return 0;
 433 }
 434
 435 #if SELFTEST
 436
 437 static char *s[] = {
 438
 439         "2^(1/2)",
 440         "2^(1/2)",
 441
 442         "2^(3/2)",
 443         "2*2^(1/2)",
 444
 445         "(-2)^(1/2)",
 446         "i*2^(1/2)",
 447
 448         "3^(4/3)",
 449         "3*3^(1/3)",
 450
 451         "3^(-4/3)",
 452         "1/(3*3^(1/3))",
 453
 454         "3^(5/3)",
 455         "3*3^(2/3)",
 456
 457         "3^(2/3)-9^(1/3)",
 458         "0",
 459
 460         "3^(10/3)",
 461         "27*3^(1/3)",
 462
 463         "3^(-10/3)",
 464         "1/(27*3^(1/3))",
 465
 466         "(1/3)^(10/3)",
 467         "1/(27*3^(1/3))",
 468
 469         "(1/3)^(-10/3)",
 470         "27*3^(1/3)",
 471
 472         "27^(2/3)",
 473         "9",
 474
 475         "27^(-2/3)",
 476         "1/9",
 477
 478         "102^(1/2)",
 479         "2^(1/2)*3^(1/2)*17^(1/2)",
 480
 481         "32^(1/3)",
 482         "2*2^(2/3)",
 483
 484         "9999^(1/2)",
 485         "3*11^(1/2)*101^(1/2)",
 486
 487         "10000^(1/3)",
 488         "10*2^(1/3)*5^(1/3)",
 489
 490         "sqrt(1000000)",
 491         "1000",
 492
 493         "sqrt(-1000000)",
 494         "1000*i",
 495
 496         "sqrt(2^60)",
 497         "1073741824",
 498
 499         // this is why we factor irrationals
 500
 501         "6^(1/3) 3^(2/3)",
 502         "3*2^(1/3)",
 503
 504         // inverse of complex numbers
 505
 506         "1/(2+3*i)",
 507         "2/13-3/13*i",
 508
 509         "1/(2+3*i)^2",
 510         "-5/169-12/169*i",
 511
 512         "(-1+3i)/(2-i)",
 513         "-1+i",
 514
 515         // other
 516
 517         "(0.0)^(0.0)",
 518         "1",
 519
 520         "(-4.0)^(1.5)",
 521         "-8*i",
 522
 523         "(-4.0)^(0.5)",
 524         "2*i",
 525
 526         "(-4.0)^(-0.5)",
 527         "-0.5*i",
 528
 529         "(-4.0)^(-1.5)",
 530         "0.125*i",
 531
 532         // more complex number cases
 533
 534         "(1+i)^2",
 535         "2*i",
 536
 537         "(1+i)^(-2)",
 538         "-1/2*i",
 539
 540         "(1+i)^(1/2)",
 541         "(-1)^(1/8)*2^(1/4)",
 542
 543         "(1+i)^(-1/2)",
 544         "-(-1)^(7/8)/(2^(1/4))",
 545
 546         "(1+i)^(0.5)",
 547         "1.09868+0.45509*i",
 548
 549         "(1+i)^(-0.5)",
 550         "0.776887-0.321797*i",
 551
 552         // test cases for simplification of polar forms, counterclockwise
 553
 554         "exp(i*pi/2)",
 555         "i",
 556
 557         "exp(i*pi)",
 558         "-1",
 559
 560         "exp(i*3*pi/2)",
 561         "-i",
 562
 563         "exp(i*2*pi)",
 564         "1",
 565
 566         "exp(i*5*pi/2)",
 567         "i",
 568
 569         "exp(i*3*pi)",
 570         "-1",
 571
 572         "exp(i*7*pi/2)",
 573         "-i",
 574
 575         "exp(i*4*pi)",
 576         "1",
 577
 578         "exp(A+i*pi/2)",
 579         "i*exp(A)",
 580
 581         "exp(A+i*pi)",
 582         "-exp(A)",
 583
 584         "exp(A+i*3*pi/2)",
 585         "-i*exp(A)",
 586
 587         "exp(A+i*2*pi)",
 588         "exp(A)",
 589
 590         "exp(A+i*5*pi/2)",
 591         "i*exp(A)",
 592
 593         "exp(A+i*3*pi)",
 594         "-exp(A)",
 595
 596         "exp(A+i*7*pi/2)",
 597         "-i*exp(A)",
 598
 599         "exp(A+i*4*pi)",
 600         "exp(A)",
 601
 602         // test cases for simplification of polar forms, clockwise
 603
 604         "exp(-i*pi/2)",
 605         "-i",
 606
 607         "exp(-i*pi)",
 608         "-1",
 609
 610         "exp(-i*3*pi/2)",
 611         "i",
 612
 613         "exp(-i*2*pi)",
 614         "1",
 615
 616         "exp(-i*5*pi/2)",
 617         "-i",
 618
 619         "exp(-i*3*pi)",
 620         "-1",
 621
 622         "exp(-i*7*pi/2)",
 623         "i",
 624
 625         "exp(-i*4*pi)",
 626         "1",
 627
 628         "exp(A-i*pi/2)",
 629         "-i*exp(A)",
 630
 631         "exp(A-i*pi)",
 632         "-exp(A)",
 633
 634         "exp(A-i*3*pi/2)",
 635         "i*exp(A)",
 636
 637         "exp(A-i*2*pi)",
 638         "exp(A)",
 639
 640         "exp(A-i*5*pi/2)",
 641         "-i*exp(A)",
 642
 643         "exp(A-i*3*pi)",
 644         "-exp(A)",
 645
 646         "exp(A-i*7*pi/2)",
 647         "i*exp(A)",
 648
 649         "exp(A-i*4*pi)",
 650         "exp(A)",
 651 };
 652
 653 void
 654 test_power(void)
 655 {
 656         test(__FILE__, s, sizeof s / sizeof (char *));
 657 }
 658
 659 #endif