From d0f686727df046c1ad0803d6eae4adbc5d3aad39 Mon Sep 17 00:00:00 2001
From: Holger Vogt <holger.vogt@uni-due.de>
Date: Sat, 12 Nov 2022 14:35:01 +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/frontend/numparam/xpressn.c | 26 +++++++++++++++++++++++++-
 1 file changed, 25 insertions(+), 1 deletion(-)

diff --git a/src/frontend/numparam/xpressn.c b/src/frontend/numparam/xpressn.c
index 3ce568223..f4732ef74 100644
--- a/src/frontend/numparam/xpressn.c
+++ b/src/frontend/numparam/xpressn.c
@@ -773,7 +773,31 @@ operate(char op, double x, double y)
         x = x / y;
         break;
     case '^':                   /* power */
-        x = pow(fabs(x), y);
+        if (newcompat.hs) {
+            if (x < 0)
+                x = pow(x, round(y));
+            else if (x == 0)
+                x = 0;
+            else
+                x = pow(x, y);
+        }
+        else if (newcompat.lt) {
+            if (x >= 0)
+                x = pow(x, y);
+            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(y), y, 10))
+                    x = pow(x, round(y));
+                else
+                    /* As per LTSPICE specification for ** */
+                    x = 0;
+            }
+        }
+        else {
+            x = pow(fabs(x), y);
+        }
         break;
     case 'A':                   /* && */
         x = ((x != 0.0) && (y != 0.0)) ? 1.0 : 0.0;