From aa7c7190dfd65b3b56b91ccd1e888dcb9498a6d4 Mon Sep 17 00:00:00 2001
From: Johannes Pfeifer <jpfeifer@gmx.de>
Date: Mon, 5 May 2025 19:04:11 +0200
Subject: [PATCH] :bug: getPowerDeriv.m: fix first-order derivative of 0^0

Evaluated to NaN instead of 0
---
 matlab/getPowerDeriv.m | 23 +++++++++++++++++++++--
 1 file changed, 21 insertions(+), 2 deletions(-)

diff --git a/matlab/getPowerDeriv.m b/matlab/getPowerDeriv.m
index 3e0d1823e5..4dee1c4172 100644
--- a/matlab/getPowerDeriv.m
+++ b/matlab/getPowerDeriv.m
@@ -30,7 +30,7 @@ function dxp=getPowerDeriv(x,p,k)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-if (abs(x) < 1e-12) && (p > 0) && (k > p) && (abs(p - round(p)) < 1e-12)
+if (abs(x) < 1e-12) && (p >= 0) && (k > p) && (abs(p - round(p)) < 1e-12)
     dxp = 0;
 else
     dxp = x^(p-k);
@@ -39,4 +39,23 @@ else
         p = p-1;
     end
 end
-end
+
+return % --*-- Unit tests --*--
+
+%@test:1
+x=getPowerDeriv(2,3,1);
+t(1)=all(abs(x-3*4)<1e-10)
+x=getPowerDeriv(0,2,2);
+t(2)=all(abs(x-2)<1e-10)
+x=getPowerDeriv(0,2,3); %special case evaluates to 0
+t(3)=all(abs(x-0)<1e-10)
+x=getPowerDeriv(1e-13,2,3-1e-13); %0 within tolerance
+t(4)=all(abs(x-0)<1e-10)
+x=getPowerDeriv(0,0,1);
+t(5)=all(abs(x-0)<1e-10)
+x=getPowerDeriv(0,0,0);
+t(6)=all(abs(x-1)<1e-10);
+x=getPowerDeriv(0,1/3,1); %derivative evaluating to Inf due to division by 0
+t(7)= isinf(x)
+T = all(t);
+%@eof:1
\ No newline at end of file
-- 
GitLab