代码拉取完成，页面将自动刷新
/* 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();
// this && n != 0x80000000 added by DDC
// as it's not always the case that 0x80000000
// is negative
if (n > 1 && ((unsigned int)n) != 0x80000000) {
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 const 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
此处可能存在不合适展示的内容，页面不予展示。您可通过相关编辑功能自查并修改。
如您确认内容无涉及 不当用语 / 纯广告导流 / 暴力 / 低俗色情 / 侵权 / 盗版 / 虚假 / 无价值内容或违法国家有关法律法规的内容，可点击提交进行申诉，我们将尽快为您处理。
马建仓 AI 助手