These routines are used on msvcrt.dll on arm, which only provides
log(). For exact powers of two, produce an exact result instead
of one that only is numerically close.

Defer to the normal implementation for anything other than the
simple cases (negative values, denormals, infinities, etc).

Signed-off-by: Martin Storsjö <[email protected]>
---
 mingw-w64-crt/math/arm-common/log2.c | 29 ++++++++++++++++++++++++++++
 1 file changed, 29 insertions(+)

diff --git a/mingw-w64-crt/math/arm-common/log2.c 
b/mingw-w64-crt/math/arm-common/log2.c
index d4fc031e9..8bf8a2412 100644
--- a/mingw-w64-crt/math/arm-common/log2.c
+++ b/mingw-w64-crt/math/arm-common/log2.c
@@ -5,14 +5,43 @@
  */
 
 #include <math.h>
+#include <stdint.h>
+
+typedef union ieee754_double_ {
+  struct __attribute__((__packed__)) {
+    uint64_t f52 : 52;
+    uint64_t exp : 11;
+    uint64_t sgn :  1;
+  };
+  double f;
+} ieee754_double;
+
+typedef union ieee754_float_ {
+  struct __attribute__((__packed__)) {
+    uint32_t f23 : 23;
+    uint32_t exp :  8;
+    uint32_t sgn :  1;
+  };
+  float f;
+} ieee754_float;
 
 double log2(double x)
 {
+    ieee754_double u = { .f = x };
+    if (u.sgn == 0 && u.f52 == 0 && u.exp > 0 && u.exp < 0x7ff) {
+        // Handle exact powers of two exactly
+        return u.exp - 1023;
+    }
     return log(x) / 0.69314718246459960938;
 }
 
 float log2f(float x)
 {
+    ieee754_float u = { .f = x };
+    if (u.sgn == 0 && u.f23 == 0 && u.exp > 0 && u.exp < 0xff) {
+        // Handle exact powers of two exactly
+        return u.exp - 127;
+    }
     return logf(x) / 0.69314718246459960938f;
 }
 
-- 
2.17.1



_______________________________________________
Mingw-w64-public mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public

Reply via email to