Add new functions for operators x**y or x^y

compatmode hs: x>0 pow(x, y), x<0 pow(x, round(y)), X=0 0
compatmode lt: x>0 pow(x, y), x<0 pow(x, y) if y is close to integer, else 0
This commit is contained in:
Holger Vogt 2022-11-12 14:35:26 +01:00
parent d0f686727d
commit f13aa89626
3 changed files with 39 additions and 2 deletions

View File

@ -123,7 +123,7 @@ static struct op {
PT_MINUS, "-", (void(*)(void)) PTminus}, {
PT_TIMES, "*", (void(*)(void)) PTtimes}, {
PT_DIVIDE, "/", (void(*)(void)) PTdivide}, {
PT_POWER, "^", (void(*)(void)) PTpower}
PT_POWER, "^", (void(*)(void)) PTpowerH}
};
#define NUM_OPS (int)NUMELEMS(ops)
@ -330,7 +330,7 @@ static INPparseNode *PTdifferentiate(INPparseNode * p, int varnum)
#define b p->right
if (b->type == PT_CONSTANT) {
arg1 = PTdifferentiate(a, varnum);
if (newcompat.lt) {
if (newcompat.hs || newcompat.lt) {
newp = mkb(PT_TIMES,
mkb(PT_TIMES,
mkcon(b->constant),

View File

@ -43,6 +43,7 @@ double PTminus(double arg1, double arg2);
double PTtimes(double arg1, double arg2);
double PTdivide(double arg1, double arg2);
double PTpower(double arg1, double arg2);
double PTpowerH(double arg1, double arg2);
double PTpwr(double arg1, double arg2);
double PTacos(double arg);
double PTacosh(double arg);

View File

@ -90,6 +90,42 @@ PTpower(double arg1, double arg2)
return res;
}
double
PTpowerH(double arg1, double arg2)
{
double res;
if (newcompat.hs) {
if (arg1 < 0)
res = pow(arg1, round(arg2));
else if (arg1 == 0){
res = 0;
}
else
{
res = pow(arg1, arg2);
}
}
else if (newcompat.lt) {
if (arg1 >= 0)
res = pow(arg1, arg2);
else {
/* If arg2 is quasi an integer, round it to have pow not fail
when arg1 is negative. Takes into account the double
representation which sometimes differs in the last digit(s). */
if (AlmostEqualUlps(nearbyint(arg2), arg2, 10))
res = pow(arg1, round(arg2));
else
/* As per LTSPICE specification for ** */
res = 0;
}
}
else
res = pow(fabs(arg1), arg2);
return res;
}
double
PTpwr(double arg1, double arg2)
{