From f13aa896265c8861d9130999baafd81a8b6f4bd4 Mon Sep 17 00:00:00 2001
From: Holger Vogt <holger.vogt@uni-due.de>
Date: Sat, 12 Nov 2022 14:35:26 +0100
Subject: [PATCH] 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

---
 src/spicelib/parser/inpptree.c |  4 ++--
 src/spicelib/parser/inpxx.h    |  1 +
 src/spicelib/parser/ptfuncs.c  | 36 ++++++++++++++++++++++++++++++++++
 3 files changed, 39 insertions(+), 2 deletions(-)

diff --git a/src/spicelib/parser/inpptree.c b/src/spicelib/parser/inpptree.c
index ca870cffb..4d725077b 100644
--- a/src/spicelib/parser/inpptree.c
+++ b/src/spicelib/parser/inpptree.c
@@ -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),
diff --git a/src/spicelib/parser/inpxx.h b/src/spicelib/parser/inpxx.h
index afbb9d162..77e4a2143 100644
--- a/src/spicelib/parser/inpxx.h
+++ b/src/spicelib/parser/inpxx.h
@@ -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);
diff --git a/src/spicelib/parser/ptfuncs.c b/src/spicelib/parser/ptfuncs.c
index 8f1dc3760..a4ed04a3c 100644
--- a/src/spicelib/parser/ptfuncs.c
+++ b/src/spicelib/parser/ptfuncs.c
@@ -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)
 {