|
|
|
@ -296,12 +296,28 @@ static INPparseNode *PTdifferentiate(INPparseNode * p, int varnum) |
|
|
|
break; |
|
|
|
|
|
|
|
case PT_POWER: |
|
|
|
/* |
|
|
|
* ^ : a^b -> |a| math^ b |
|
|
|
* |
|
|
|
* D(pow(a,b)) |
|
|
|
* = D(exp(b*log(abs(a)))) |
|
|
|
* = exp(b*log(abs(a))) * D(b*log(abs(a))) |
|
|
|
* = pow(a,b) * (D(b)*log(abs(a)) + b*D(abs(a))/abs(a)) |
|
|
|
* = pow(a,b) * (D(b)*log(abs(a)) + b*sgn(a)*D(a)/abs(a)) |
|
|
|
* = pow(a,b) * (D(b)*log(abs(a)) + b*D(a)/a) |
|
|
|
* |
|
|
|
* when D(b) == 0, then |
|
|
|
* |
|
|
|
* D(pow(a,b)) |
|
|
|
* = pow(a,b) * (D(b)*log(abs(a)) + b*D(a)/a) |
|
|
|
* = pow(a,b) * b * D(a)/a |
|
|
|
* = pow(a,b) * b * D(a)/(signum(a) * abs(a)) |
|
|
|
* = pow(a, b-1) * b * D(a) / signum(a) |
|
|
|
* = pwr(a, b-1) * b * D(a) |
|
|
|
*/ |
|
|
|
#define a p->left |
|
|
|
#define b p->right |
|
|
|
if (b->type == PT_CONSTANT) { |
|
|
|
/* |
|
|
|
* D(a^C) = C * a^(C-1) * D(a) |
|
|
|
*/ |
|
|
|
arg1 = PTdifferentiate(a, varnum); |
|
|
|
|
|
|
|
newp = mkb(PT_TIMES, mkb(PT_TIMES, |
|
|
|
@ -310,11 +326,6 @@ static INPparseNode *PTdifferentiate(INPparseNode * p, int varnum) |
|
|
|
mkcon(b->constant - 1))), |
|
|
|
arg1); |
|
|
|
} else { |
|
|
|
/* |
|
|
|
* D(a^b) = D(exp(b*log(a))) |
|
|
|
* = exp(b*log(a)) * D(b*log(a)) |
|
|
|
* = exp(b*log(a)) * (D(b)*log(a) + b*D(a)/a) |
|
|
|
*/ |
|
|
|
arg1 = PTdifferentiate(a, varnum); |
|
|
|
arg2 = PTdifferentiate(b, varnum); |
|
|
|
newp = mkb(PT_TIMES, mkf(PTF_EXP, mkb(PT_TIMES, |
|
|
|
@ -551,6 +562,25 @@ static INPparseNode *PTdifferentiate(INPparseNode * p, int varnum) |
|
|
|
break; |
|
|
|
|
|
|
|
case PTF_POW: |
|
|
|
/* |
|
|
|
* pow : pow(a,b) -> |a| math^ b |
|
|
|
* |
|
|
|
* D(pow(a,b)) |
|
|
|
* = D(exp(b*log(abs(a)))) |
|
|
|
* = exp(b*log(abs(a))) * D(b*log(abs(a))) |
|
|
|
* = pow(a,b) * (D(b)*log(abs(a)) + b*D(abs(a))/abs(a)) |
|
|
|
* = pow(a,b) * (D(b)*log(abs(a)) + b*sgn(a)*D(a)/abs(a)) |
|
|
|
* = pow(a,b) * (D(b)*log(abs(a)) + b*D(a)/a) |
|
|
|
* |
|
|
|
* when D(b) == 0, then |
|
|
|
* |
|
|
|
* D(pow(a,b)) |
|
|
|
* = pow(a,b) * (D(b)*log(abs(a)) + b*D(a)/a) |
|
|
|
* = pow(a,b) * b * D(a)/a |
|
|
|
* = pow(a,b) * b * D(a)/(signum(a) * abs(a)) |
|
|
|
* = pow(a, b-1) * b * D(a) / signum(a) |
|
|
|
* = pwr(a, b-1) * b * D(a) |
|
|
|
*/ |
|
|
|
{ |
|
|
|
/* |
|
|
|
pow(a,b) |
|
|
|
@ -560,9 +590,6 @@ static INPparseNode *PTdifferentiate(INPparseNode * p, int varnum) |
|
|
|
#define b p->left->right |
|
|
|
|
|
|
|
if (b->type == PT_CONSTANT) { |
|
|
|
/* |
|
|
|
* D(f^C) = C * f^(C-1) * D(f) |
|
|
|
*/ |
|
|
|
arg1 = PTdifferentiate(a, varnum); |
|
|
|
|
|
|
|
newp = mkb(PT_TIMES, mkb(PT_TIMES, |
|
|
|
@ -576,11 +603,6 @@ static INPparseNode *PTdifferentiate(INPparseNode * p, int varnum) |
|
|
|
printf("\n"); |
|
|
|
#endif |
|
|
|
} else { |
|
|
|
/* |
|
|
|
* D(f^g) = D(exp(g*log(f))) |
|
|
|
* = exp(g*log(f)) * D(g*log(f)) |
|
|
|
* = exp(g*log(f)) * (D(g)*log(f) + g*D(f)/f) |
|
|
|
*/ |
|
|
|
arg1 = PTdifferentiate(a, varnum); |
|
|
|
arg2 = PTdifferentiate(b, varnum); |
|
|
|
newp = mkb(PT_TIMES, mkf(PTF_EXP, mkb(PT_TIMES, |
|
|
|
@ -599,6 +621,26 @@ static INPparseNode *PTdifferentiate(INPparseNode * p, int varnum) |
|
|
|
break; |
|
|
|
|
|
|
|
case PTF_PWR: |
|
|
|
/* |
|
|
|
* pwr : pwr(a,b) -> signum(a) * (|a| math^ b) |
|
|
|
* -> signum(a) * pow(a, b) |
|
|
|
* |
|
|
|
* Note: |
|
|
|
* D(pow(a,b)) = pow(a,b) * (D(b)*log(abs(a)) + b*D(a)/a) |
|
|
|
* |
|
|
|
* D(pwr(a,b)) |
|
|
|
* = D(signum(a) * pow(a,b)) |
|
|
|
* = D(signum(a)) * pow(a,b) + signum(a) * D(pow(a,b)) |
|
|
|
* = 0 + signum(a) * pow(a,b) * (D(b)*log(abs(a)) + b*D(a)/a) |
|
|
|
* = pwr(a,b) * (D(b)*log(abs(a)) + b*D(a)/a) |
|
|
|
* |
|
|
|
* with D(b) == 0 |
|
|
|
* |
|
|
|
* D(pwr(a,b)) |
|
|
|
* = pwr(a,b) * b * D(a)/a |
|
|
|
* = signum(a) * pow(a,b) * b * D(a)/(signum(a) * abs(a)) |
|
|
|
* = pow(a, b-1) * b * D(a) |
|
|
|
*/ |
|
|
|
{ |
|
|
|
/* |
|
|
|
pwr(a,b) |
|
|
|
@ -607,18 +649,6 @@ static INPparseNode *PTdifferentiate(INPparseNode * p, int varnum) |
|
|
|
#define a p->left->left |
|
|
|
#define b p->left->right |
|
|
|
if (b->type == PT_CONSTANT) { |
|
|
|
/* b is a constant |
|
|
|
* |
|
|
|
* f(a,b) = signum(a) * abs(a)^b |
|
|
|
* = signum(a) * exp(b*log(abs(a))) |
|
|
|
* D(f) = signum(a) * D(exp(b*log(abs(a)))) |
|
|
|
* = signum(a) * exp(b*log(abs(a))) * D(b*log(abs(a))) |
|
|
|
* = signum(a) * abs(a)^b * D(b*log(abs(a))) |
|
|
|
* = signum(a) * abs(a)^b * b * 1/abs(a) * D(abs(a)) |
|
|
|
* = signum(a) * abs(a)^(b-1) * b * D(abs(a)) |
|
|
|
* = signum(a) * abs(a)^(b-1) * b * signum(a) * D(a) |
|
|
|
* = abs(a)^(b-1) * b * D(a) |
|
|
|
*/ |
|
|
|
arg1 = PTdifferentiate(a, varnum); |
|
|
|
|
|
|
|
newp = mkb(PT_TIMES, |
|
|
|
@ -634,19 +664,6 @@ static INPparseNode *PTdifferentiate(INPparseNode * p, int varnum) |
|
|
|
printf("\n"); |
|
|
|
#endif |
|
|
|
} else { |
|
|
|
/* b is a function |
|
|
|
* |
|
|
|
* f(a,b) = signum(a) * abs(a)^b |
|
|
|
* = signum(a) * exp(b*log(abs(a))) |
|
|
|
* D(f) = signum(a) * D(exp(b*log(abs(a)))) |
|
|
|
* = signum(a) * exp(b*log(abs(a))) * D(b*log(abs(a))) |
|
|
|
* = signum(a) * exp(b*log(abs(a))) * (D(b) * log(abs(a)) + b * D(log(abs(a)))) |
|
|
|
* = signum(a) * exp(b*log(abs(a))) * (D(b) * log(abs(a)) + b * 1/abs(a) * D(abs(a))) |
|
|
|
* = signum(a) * exp(b*log(abs(a))) * (D(b) * log(abs(a)) + b * 1/abs(a) * signum(a)*D(a)) |
|
|
|
* = signum(a) * exp(b*log(abs(a))) * (D(b) * log(abs(a)) + b/a*D(a)) |
|
|
|
* = signum(a) * exp(b*log(abs(a))) * D(b) * log(abs(a) + signum(a) * exp(b*log(abs(a))) / a * b * D(a) |
|
|
|
* = signum(a) * exp(b*log(abs(a))) * D(b) * log(abs(a) + abs(a)^(b-1) * b * D(a) |
|
|
|
*/ |
|
|
|
arg1 = PTdifferentiate(a, varnum); |
|
|
|
arg2 = PTdifferentiate(b, varnum); |
|
|
|
newp = mkb(PT_PLUS, |
|
|
|
|