Changeset 516e780 in mainline for uspace/lib/math/generic/round.c

Timestamp:

2018-08-31T11:55:41Z (6 years ago)

Author:

GitHub <noreply@…>

Branches:

lfn, master, serial, ticket/834-toolchain-update, topic/msim-upgrade, topic/simplify-dev-export

Children:

fa86fff

Parents:

7f7d642

git-author:

Jiří Zárevúcky <zarevucky.jiri@…> (2018-08-31 11:55:41)

git-committer:

GitHub <noreply@…> (2018-08-31 11:55:41)

Message:

Strip down libmath. (#45)

libmath is mostly unused (except for trunc(), sin() and cos()), and most functions in it are either very imprecise or downright broken. Additionally, it is implemented in manner that conflicts with C standard. Instead of trying to fix all the shortcomings while maintaining unused functionality, I'm opting to simply remove most of it and only keep the parts that are currently necessary.

Later readdition of the removed functions is possible, but there needs to be a reliable way to evaluate their quality first.

File:

• uspace/lib/math/generic/round.c (moved) (moved from uspace/lib/math/generic/pow.c ) (1 diff)

uspace/lib/math/generic/round.c

@@ -35 +35 @@
 
 #include <math.h>
-#include <pow.h>
+#include <float.h>
+#include <stdint.h>
 
-/** Single precision power
- *
- * Compute power value.
- *
- * @param x Base
- * @param y Exponent
- *
- * @return Cosine value.
- *
+/**
+ * Rounds its argument to the nearest integer value in floating-point format,
+ * rounding halfway cases away from zero, regardless of the current rounding
+ * direction.
  */
-float32_t float32_pow(float32_t x, float32_t y)
+float roundf(float val)
 {
-        /* x^y = (e ^ log(x))^y = e ^ (log(x) * y) */
-        return exp_f32(log_f32(x) * y);
+        const int exp_bias = FLT_MAX_EXP - 1;
+        const int mant_bits = FLT_MANT_DIG - 1;
+
+        union {
+                float f;
+                uint32_t i;
+        } u = { .f = fabsf(val) };
+
+        int exp = (u.i >> mant_bits) - exp_bias;
+
+        /* If value is less than 0.5, return zero with appropriate sign. */
+        if (exp < -1)
+                return copysignf(0.0f, val);
+
+        /* If exponent is exactly mant_bits, adding 0.5 could change the result. */
+        if (exp >= mant_bits)
+                return val;
+
+        /* Use trunc with adjusted value to do the rounding. */
+        return copysignf(truncf(u.f + 0.5), val);
 }
 
-/** Double precision power
- *
- * Compute power value.
- *
- * @param x Base
- * @param y Exponent
- *
- * @return Cosine value.
- *
+/**
+ * Rounds its argument to the nearest integer value in floating-point format,
+ * rounding halfway cases away from zero, regardless of the current rounding
+ * direction.
  */
-float64_t float64_pow(float64_t x, float64_t y)
+double round(double val)
 {
-        /* x^y = (e ^ log(x))^y = e ^ (log(x) * y) */
-        return exp_f64(log_f64(x) * y);
+        const int exp_bias = DBL_MAX_EXP - 1;
+        const int mant_bits = DBL_MANT_DIG - 1;
+
+        union {
+                double f;
+                uint64_t i;
+        } u = { .f = fabs(val) };
+
+        int exp = ((int)(u.i >> mant_bits)) - exp_bias;
+
+        /* If value is less than 0.5, return zero with appropriate sign. */
+        if (exp < -1)
+                return copysign(0.0, val);
+
+        /* If exponent is exactly mant_bits, adding 0.5 could change the result. */
+        if (exp >= mant_bits)
+                return val;
+
+        /* Use trunc with adjusted value to do the rounding. */
+        return copysign(trunc(u.f + 0.5), val);
 }