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]>
---
Added a cast to int before returning the adjusted exponent as a
float.
---
 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..42ae3eeeb 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 (int)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 (int)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