[libm] Rough approximations of expf, expm1f, logf, log1pf

esbmc · lucasccordeiro · Jan 19, 2024 · Nov 16, 2023 · Nov 16, 2023 · Nov 16, 2023
commit 02472d0a35323995e65a0f541496270f7007f3d6
diff --git a/src/c2goto/library/libm/expf.c b/src/c2goto/library/libm/expf.c
@@ -0,0 +1,42 @@
+#include <math.h>
+
+float expf(float x)
+{
+__ESBMC_HIDE:;
+    /**
+        Highly imprecise implementation of expf(x) on floating-point operands
+        which satisfies the following properties:
+        expf(x) >= 0
+        expf(x) >= 1 + x
+        expf(-Inf) = 0
+        expf(0) = 1
+        expf(+Inf) = Inf
+        dexpf(x=0)/dx = 1
+        logf(expf(x)) = x
+     */
+
+    // handle NaN inputs
+    if(isnan(x))
+        return x;
+
+    // handle infinities (x*x will overflow before x)
+    float xx = x * x;
+    if(isinf(xx)) {
+        if(x < 0.0f)
+            return 0.0f; // expf(-Inf) = 0
+        else
+            return 1.0f / 0.0f; // expf(+Inf) = Inf
+    }
+
+    // inverse function of the logf approximation
+    // i.e. inverse((x - 1) / sqrt(x))
+    // it is numerically unstable for x < 0
+    // since it relies on x - sqrt(x^2 + 4) cancelling out
+    return 0.5 * x * (x + sqrtf(xx + 4.0f)) + 1.0f;
+}
+
+float expm1f(float x)
+{
+__ESBMC_HIDE:;
+    return expf(x) - 1.0f;
+}
diff --git a/src/c2goto/library/libm/logf.c b/src/c2goto/library/libm/logf.c
@@ -0,0 +1,33 @@
+#include <math.h>
+
+float logf(float x)
+{
+__ESBMC_HIDE:;
+    /**
+        Highly imprecise implementation of logf(x) on floating-point operands
+        which satisfies the following properties:
+        logf(x) <= x - 1
+        logf(0) = -Inf
+        logf(1) = 0
+        logf(Inf) = Inf
+        dlogf(x=1)/dx = 1
+        expf(logf(x)) = x
+     */
+
+    // handle negative inputs
+    if(x < 0.0f)
+        return 0.0f / 0.0f; // return NaN
+
+    // handle positive infinity
+    if(isinf(x))
+        return x;
+
+    // ad-hoc fractional Pade' approximant
+    return (x - 1.0f) / sqrtf(x);
+}
+
+float log1pf(float x)
+{
+__ESBMC_HIDE:;
+    return logf(x + 1.0f);
+}