660 lines
8.3 KiB
C++
660 lines
8.3 KiB
C++
/* Power function
|
|
|
|
Input: push Base
|
|
|
|
push Exponent
|
|
|
|
Output: Result on stack
|
|
*/
|
|
|
|
#include "stdafx.h"
|
|
#include "defs.h"
|
|
|
|
void
|
|
eval_power(void)
|
|
{
|
|
push(cadr(p1));
|
|
eval();
|
|
push(caddr(p1));
|
|
eval();
|
|
power();
|
|
}
|
|
|
|
void
|
|
power(void)
|
|
{
|
|
save();
|
|
yypower();
|
|
restore();
|
|
}
|
|
|
|
void
|
|
yypower(void)
|
|
{
|
|
int n;
|
|
|
|
p2 = pop();
|
|
p1 = pop();
|
|
|
|
// both base and exponent are rational numbers?
|
|
|
|
if (isrational(p1) && isrational(p2)) {
|
|
push(p1);
|
|
push(p2);
|
|
qpow();
|
|
return;
|
|
}
|
|
|
|
// both base and exponent are either rational or double?
|
|
|
|
if (isnum(p1) && isnum(p2)) {
|
|
push(p1);
|
|
push(p2);
|
|
dpow();
|
|
return;
|
|
}
|
|
|
|
if (istensor(p1)) {
|
|
power_tensor();
|
|
return;
|
|
}
|
|
|
|
if (p1 == symbol(E) && car(p2) == symbol(LOG)) {
|
|
push(cadr(p2));
|
|
return;
|
|
}
|
|
|
|
if (p1 == symbol(E) && isdouble(p2)) {
|
|
push_double(exp(p2->u.d));
|
|
return;
|
|
}
|
|
|
|
// 1 ^ a -> 1
|
|
|
|
// a ^ 0 -> 1
|
|
|
|
if (equal(p1, one) || iszero(p2)) {
|
|
push(one);
|
|
return;
|
|
}
|
|
|
|
// a ^ 1 -> a
|
|
|
|
if (equal(p2, one)) {
|
|
push(p1);
|
|
return;
|
|
}
|
|
|
|
// (a * b) ^ c -> (a ^ c) * (b ^ c)
|
|
|
|
if (car(p1) == symbol(MULTIPLY)) {
|
|
p1 = cdr(p1);
|
|
push(car(p1));
|
|
push(p2);
|
|
power();
|
|
p1 = cdr(p1);
|
|
while (iscons(p1)) {
|
|
push(car(p1));
|
|
push(p2);
|
|
power();
|
|
multiply();
|
|
p1 = cdr(p1);
|
|
}
|
|
return;
|
|
}
|
|
|
|
// (a ^ b) ^ c -> a ^ (b * c)
|
|
|
|
if (car(p1) == symbol(POWER)) {
|
|
push(cadr(p1));
|
|
push(caddr(p1));
|
|
push(p2);
|
|
multiply();
|
|
power();
|
|
return;
|
|
}
|
|
|
|
// (a + b) ^ n -> (a + b) * (a + b) ...
|
|
|
|
if (expanding && isadd(p1) && isnum(p2)) {
|
|
push(p2);
|
|
n = pop_integer();
|
|
if (n > 1) {
|
|
power_sum(n);
|
|
return;
|
|
}
|
|
}
|
|
|
|
// sin(x) ^ 2n -> (1 - cos(x) ^ 2) ^ n
|
|
|
|
if (trigmode == 1 && car(p1) == symbol(SIN) && iseveninteger(p2)) {
|
|
push_integer(1);
|
|
push(cadr(p1));
|
|
cosine();
|
|
push_integer(2);
|
|
power();
|
|
subtract();
|
|
push(p2);
|
|
push_rational(1, 2);
|
|
multiply();
|
|
power();
|
|
return;
|
|
}
|
|
|
|
// cos(x) ^ 2n -> (1 - sin(x) ^ 2) ^ n
|
|
|
|
if (trigmode == 2 && car(p1) == symbol(COS) && iseveninteger(p2)) {
|
|
push_integer(1);
|
|
push(cadr(p1));
|
|
sine();
|
|
push_integer(2);
|
|
power();
|
|
subtract();
|
|
push(p2);
|
|
push_rational(1, 2);
|
|
multiply();
|
|
power();
|
|
return;
|
|
}
|
|
|
|
// complex number? (just number, not expression)
|
|
|
|
if (iscomplexnumber(p1)) {
|
|
|
|
// integer power?
|
|
|
|
// n will be negative here, positive n already handled
|
|
|
|
if (isinteger(p2)) {
|
|
|
|
// / \ n
|
|
// -n | a - ib |
|
|
// (a + ib) = | -------- |
|
|
// | 2 2 |
|
|
// \ a + b /
|
|
|
|
push(p1);
|
|
conjugate();
|
|
p3 = pop();
|
|
push(p3);
|
|
push(p3);
|
|
push(p1);
|
|
multiply();
|
|
divide();
|
|
push(p2);
|
|
negate();
|
|
power();
|
|
return;
|
|
}
|
|
|
|
// noninteger or floating power?
|
|
|
|
if (isnum(p2)) {
|
|
|
|
#if 1 // use polar form
|
|
push(p1);
|
|
mag();
|
|
push(p2);
|
|
power();
|
|
push_integer(-1);
|
|
push(p1);
|
|
arg();
|
|
push(p2);
|
|
multiply();
|
|
push(symbol(PI));
|
|
divide();
|
|
power();
|
|
multiply();
|
|
|
|
#else // use exponential form
|
|
push(p1);
|
|
mag();
|
|
push(p2);
|
|
power();
|
|
push(symbol(E));
|
|
push(p1);
|
|
arg();
|
|
push(p2);
|
|
multiply();
|
|
push(imaginaryunit);
|
|
multiply();
|
|
power();
|
|
multiply();
|
|
#endif
|
|
return;
|
|
}
|
|
}
|
|
|
|
if (simplify_polar())
|
|
return;
|
|
|
|
push_symbol(POWER);
|
|
push(p1);
|
|
push(p2);
|
|
list(3);
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
//
|
|
// Compute the power of a sum
|
|
//
|
|
// Input: p1 sum
|
|
//
|
|
// n exponent
|
|
//
|
|
// Output: Result on stack
|
|
//
|
|
// Note:
|
|
//
|
|
// Uses the multinomial series (see Math World)
|
|
//
|
|
// n n! n1 n2 nk
|
|
// (a1 + a2 + ... + ak) = sum (--------------- a1 a2 ... ak )
|
|
// n1! n2! ... nk!
|
|
//
|
|
// The sum is over all n1 ... nk such that n1 + n2 + ... + nk = n.
|
|
//
|
|
//-----------------------------------------------------------------------------
|
|
|
|
// first index is the term number 0..k-1, second index is the exponent 0..n
|
|
|
|
#define A(i, j) frame[(i) * (n + 1) + (j)]
|
|
|
|
void
|
|
power_sum(int n)
|
|
{
|
|
int *a, i, j, k;
|
|
|
|
// number of terms in the sum
|
|
|
|
k = length(p1) - 1;
|
|
|
|
// local frame
|
|
|
|
push_frame(k * (n + 1));
|
|
|
|
// array of powers
|
|
|
|
p1 = cdr(p1);
|
|
for (i = 0; i < k; i++) {
|
|
for (j = 0; j <= n; j++) {
|
|
push(car(p1));
|
|
push_integer(j);
|
|
power();
|
|
A(i, j) = pop();
|
|
}
|
|
p1 = cdr(p1);
|
|
}
|
|
|
|
push_integer(n);
|
|
factorial();
|
|
p1 = pop();
|
|
|
|
a = (int *) malloc(k * sizeof (int));
|
|
|
|
if (a == NULL)
|
|
stop("malloc failure");
|
|
|
|
push(zero);
|
|
|
|
multinomial_sum(k, n, a, 0, n);
|
|
|
|
free(a);
|
|
|
|
pop_frame(k * (n + 1));
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
//
|
|
// Compute multinomial sum
|
|
//
|
|
// Input: k number of factors
|
|
//
|
|
// n overall exponent
|
|
//
|
|
// a partition array
|
|
//
|
|
// i partition array index
|
|
//
|
|
// m partition remainder
|
|
//
|
|
// p1 n!
|
|
//
|
|
// A factor array
|
|
//
|
|
// Output: Result on stack
|
|
//
|
|
// Note:
|
|
//
|
|
// Uses recursive descent to fill the partition array.
|
|
//
|
|
//-----------------------------------------------------------------------------
|
|
|
|
void
|
|
multinomial_sum(int k, int n, int *a, int i, int m)
|
|
{
|
|
int j;
|
|
|
|
if (i < k - 1) {
|
|
for (j = 0; j <= m; j++) {
|
|
a[i] = j;
|
|
multinomial_sum(k, n, a, i + 1, m - j);
|
|
}
|
|
return;
|
|
}
|
|
|
|
a[i] = m;
|
|
|
|
// coefficient
|
|
|
|
push(p1);
|
|
|
|
for (j = 0; j < k; j++) {
|
|
push_integer(a[j]);
|
|
factorial();
|
|
divide();
|
|
}
|
|
|
|
// factors
|
|
|
|
for (j = 0; j < k; j++) {
|
|
push(A(j, a[j]));
|
|
multiply();
|
|
}
|
|
|
|
add();
|
|
}
|
|
|
|
// exp(n/2 i pi) ?
|
|
|
|
// p2 is the exponent expression
|
|
|
|
// clobbers p3
|
|
|
|
int
|
|
simplify_polar(void)
|
|
{
|
|
int n;
|
|
|
|
n = isquarterturn(p2);
|
|
|
|
switch(n) {
|
|
case 0:
|
|
break;
|
|
case 1:
|
|
push_integer(1);
|
|
return 1;
|
|
case 2:
|
|
push_integer(-1);
|
|
return 1;
|
|
case 3:
|
|
push(imaginaryunit);
|
|
return 1;
|
|
case 4:
|
|
push(imaginaryunit);
|
|
negate();
|
|
return 1;
|
|
}
|
|
|
|
if (car(p2) == symbol(ADD)) {
|
|
p3 = cdr(p2);
|
|
while (iscons(p3)) {
|
|
n = isquarterturn(car(p3));
|
|
if (n)
|
|
break;
|
|
p3 = cdr(p3);
|
|
}
|
|
switch (n) {
|
|
case 0:
|
|
return 0;
|
|
case 1:
|
|
push_integer(1);
|
|
break;
|
|
case 2:
|
|
push_integer(-1);
|
|
break;
|
|
case 3:
|
|
push(imaginaryunit);
|
|
break;
|
|
case 4:
|
|
push(imaginaryunit);
|
|
negate();
|
|
break;
|
|
}
|
|
push(p2);
|
|
push(car(p3));
|
|
subtract();
|
|
exponential();
|
|
multiply();
|
|
return 1;
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
#if SELFTEST
|
|
|
|
static char *s[] = {
|
|
|
|
"2^(1/2)",
|
|
"2^(1/2)",
|
|
|
|
"2^(3/2)",
|
|
"2*2^(1/2)",
|
|
|
|
"(-2)^(1/2)",
|
|
"i*2^(1/2)",
|
|
|
|
"3^(4/3)",
|
|
"3*3^(1/3)",
|
|
|
|
"3^(-4/3)",
|
|
"1/(3*3^(1/3))",
|
|
|
|
"3^(5/3)",
|
|
"3*3^(2/3)",
|
|
|
|
"3^(2/3)-9^(1/3)",
|
|
"0",
|
|
|
|
"3^(10/3)",
|
|
"27*3^(1/3)",
|
|
|
|
"3^(-10/3)",
|
|
"1/(27*3^(1/3))",
|
|
|
|
"(1/3)^(10/3)",
|
|
"1/(27*3^(1/3))",
|
|
|
|
"(1/3)^(-10/3)",
|
|
"27*3^(1/3)",
|
|
|
|
"27^(2/3)",
|
|
"9",
|
|
|
|
"27^(-2/3)",
|
|
"1/9",
|
|
|
|
"102^(1/2)",
|
|
"2^(1/2)*3^(1/2)*17^(1/2)",
|
|
|
|
"32^(1/3)",
|
|
"2*2^(2/3)",
|
|
|
|
"9999^(1/2)",
|
|
"3*11^(1/2)*101^(1/2)",
|
|
|
|
"10000^(1/3)",
|
|
"10*2^(1/3)*5^(1/3)",
|
|
|
|
"sqrt(1000000)",
|
|
"1000",
|
|
|
|
"sqrt(-1000000)",
|
|
"1000*i",
|
|
|
|
"sqrt(2^60)",
|
|
"1073741824",
|
|
|
|
// this is why we factor irrationals
|
|
|
|
"6^(1/3) 3^(2/3)",
|
|
"3*2^(1/3)",
|
|
|
|
// inverse of complex numbers
|
|
|
|
"1/(2+3*i)",
|
|
"2/13-3/13*i",
|
|
|
|
"1/(2+3*i)^2",
|
|
"-5/169-12/169*i",
|
|
|
|
"(-1+3i)/(2-i)",
|
|
"-1+i",
|
|
|
|
// other
|
|
|
|
"(0.0)^(0.0)",
|
|
"1",
|
|
|
|
"(-4.0)^(1.5)",
|
|
"-8*i",
|
|
|
|
"(-4.0)^(0.5)",
|
|
"2*i",
|
|
|
|
"(-4.0)^(-0.5)",
|
|
"-0.5*i",
|
|
|
|
"(-4.0)^(-1.5)",
|
|
"0.125*i",
|
|
|
|
// more complex number cases
|
|
|
|
"(1+i)^2",
|
|
"2*i",
|
|
|
|
"(1+i)^(-2)",
|
|
"-1/2*i",
|
|
|
|
"(1+i)^(1/2)",
|
|
"(-1)^(1/8)*2^(1/4)",
|
|
|
|
"(1+i)^(-1/2)",
|
|
"-(-1)^(7/8)/(2^(1/4))",
|
|
|
|
"(1+i)^(0.5)",
|
|
"1.09868+0.45509*i",
|
|
|
|
"(1+i)^(-0.5)",
|
|
"0.776887-0.321797*i",
|
|
|
|
// test cases for simplification of polar forms, counterclockwise
|
|
|
|
"exp(i*pi/2)",
|
|
"i",
|
|
|
|
"exp(i*pi)",
|
|
"-1",
|
|
|
|
"exp(i*3*pi/2)",
|
|
"-i",
|
|
|
|
"exp(i*2*pi)",
|
|
"1",
|
|
|
|
"exp(i*5*pi/2)",
|
|
"i",
|
|
|
|
"exp(i*3*pi)",
|
|
"-1",
|
|
|
|
"exp(i*7*pi/2)",
|
|
"-i",
|
|
|
|
"exp(i*4*pi)",
|
|
"1",
|
|
|
|
"exp(A+i*pi/2)",
|
|
"i*exp(A)",
|
|
|
|
"exp(A+i*pi)",
|
|
"-exp(A)",
|
|
|
|
"exp(A+i*3*pi/2)",
|
|
"-i*exp(A)",
|
|
|
|
"exp(A+i*2*pi)",
|
|
"exp(A)",
|
|
|
|
"exp(A+i*5*pi/2)",
|
|
"i*exp(A)",
|
|
|
|
"exp(A+i*3*pi)",
|
|
"-exp(A)",
|
|
|
|
"exp(A+i*7*pi/2)",
|
|
"-i*exp(A)",
|
|
|
|
"exp(A+i*4*pi)",
|
|
"exp(A)",
|
|
|
|
// test cases for simplification of polar forms, clockwise
|
|
|
|
"exp(-i*pi/2)",
|
|
"-i",
|
|
|
|
"exp(-i*pi)",
|
|
"-1",
|
|
|
|
"exp(-i*3*pi/2)",
|
|
"i",
|
|
|
|
"exp(-i*2*pi)",
|
|
"1",
|
|
|
|
"exp(-i*5*pi/2)",
|
|
"-i",
|
|
|
|
"exp(-i*3*pi)",
|
|
"-1",
|
|
|
|
"exp(-i*7*pi/2)",
|
|
"i",
|
|
|
|
"exp(-i*4*pi)",
|
|
"1",
|
|
|
|
"exp(A-i*pi/2)",
|
|
"-i*exp(A)",
|
|
|
|
"exp(A-i*pi)",
|
|
"-exp(A)",
|
|
|
|
"exp(A-i*3*pi/2)",
|
|
"i*exp(A)",
|
|
|
|
"exp(A-i*2*pi)",
|
|
"exp(A)",
|
|
|
|
"exp(A-i*5*pi/2)",
|
|
"-i*exp(A)",
|
|
|
|
"exp(A-i*3*pi)",
|
|
"-exp(A)",
|
|
|
|
"exp(A-i*7*pi/2)",
|
|
"i*exp(A)",
|
|
|
|
"exp(A-i*4*pi)",
|
|
"exp(A)",
|
|
};
|
|
|
|
void
|
|
test_power(void)
|
|
{
|
|
test(__FILE__, s, sizeof s / sizeof (char *));
|
|
}
|
|
|
|
#endif
|