1

libs/libc/math/Make.defs

@@ -27,7 +27,7 @@ CSRCS += lib_coshf.c  lib_expf.c lib_fabsf.c lib_fmodf.c lib_frexpf.c
 CSRCS += lib_ldexpf.c lib_logf.c lib_log10f.c lib_log2f.c lib_modff.c
 CSRCS += lib_powf.c lib_sinf.c lib_sinhf.c lib_sqrtf.c lib_tanf.c
 CSRCS += lib_tanhf.c lib_asinhf.c lib_acoshf.c lib_atanhf.c lib_erff.c
-CSRCS += lib_copysignf.c
+CSRCS += lib_copysignf.c lib_scalbnf.c lib_scalbn.c lib_scalbnl.c
 
 CSRCS += lib_acos.c lib_asin.c lib_atan.c lib_atan2.c lib_cos.c
 CSRCS += lib_cosh.c lib_exp.c lib_fabs.c lib_fmod.c lib_frexp.c

libs/libc/math/lib_scalbn.c

@@ -0,0 +1,88 @@
+/****************************************************************************
+ * libs/libc/math/lib_scalbn.c
+ * get a double number of x*2^n
+ *
+ * This file is copy from musl libc
+ * musl is an implementation of the C standard library built on top of the
+ * Linux system call API, including interfaces defined in the base language
+ * standard, POSIX, and widely agreed-upon extensions.
+ * musl is lightweight, fast, simple, free, and strives to be correct in
+ * the sense of standards-conformance and safety.
+ *
+ *
+ ****************************************************************************/
+
+/****************************************************************************
+ * Included Files
+ ****************************************************************************/
+
+#include <math.h>
+#include <stdint.h>
+
+/****************************************************************************
+ * Pre-processor definitions
+ ****************************************************************************/
+
+#ifdef CONFIG_HAVE_DOUBLE
+#define DOUBLE_MIN (0x1p-1022)
+
+/****************************************************************************
+ * Public Functions
+ ****************************************************************************/
+
+/****************************************************************************
+ * Name: scalbn
+ *
+ * Description:
+ *   get a double number of x*2^n
+ *
+ ****************************************************************************/
+
+double scalbn(double x, int n)
+{
+  union
+    {
+      double f;
+      uint64_t i;
+    }u;
+
+  double_t y = x;
+
+  if (n > 1023)
+    {
+      y *= 0x1p1023;
+      n -= 1023;
+      if (n > 1023)
+        {
+          y *= 0x1p1023;
+          n -= 1023;
+          if (n > 1023)
+            {
+              n = 1023;
+            }
+        }
+    }
+  else if (n < -1022)
+    {
+      /* make sure final n < -53 to avoid double
+       * rounding in the subnormal range
+       */
+
+      y *= DOUBLE_MIN * 0x1p53;
+      n += 1022 - 53;
+      if (n < -1022)
+        {
+          y *= DOUBLE_MIN * 0x1p53;
+          n += 1022 - 53;
+          if (n < -1022)
+            {
+              n = -1022;
+            }
+        }
+    }
+
+  u.i = (uint64_t)(0x3ff + n) << 52;
+  x = y * u.f;
+  return x;
+}
+#endif /* CONFIG_HAVE_DOUBLE */

libs/libc/math/lib_scalbnf.c

@@ -0,0 +1,83 @@
+/****************************************************************************
+ * libs/libc/math/lib_scalbnf.c
+ * get a float number of x*2^n
+ *
+ * This file is copy from musl libc
+ * musl is an implementation of the C standard library built on top of the
+ * Linux system call API, including interfaces defined in the base language
+ * standard, POSIX, and widely agreed-upon extensions.
+ * musl is lightweight, fast, simple, free, and strives to be correct in
+ * the sense of standards-conformance and safety.
+ *
+ *
+ ****************************************************************************/
+
+/****************************************************************************
+ * Included Files
+ ****************************************************************************/
+
+#include <math.h>
+#include <stdint.h>
+
+/****************************************************************************
+ * Pre-processor definitions
+ ****************************************************************************/
+
+#define FLT_MIN (0x1p-126f)
+
+/****************************************************************************
+ * Public Functions
+ ****************************************************************************/
+
+/****************************************************************************
+ * Name: scalbnf
+ *
+ * Description:
+ *   get a float number of x*2^n
+ *
+ ****************************************************************************/
+
+float scalbnf(float x, int n)
+{
+  union
+    {
+      float f;
+      uint32_t i;
+    }u;
+
+  float_t y = x;
+
+  if (n > 127)
+    {
+      y *= 0x1p127f;
+      n -= 127;
+      if (n > 127)
+        {
+          y *= 0x1p127f;
+          n -= 127;
+          if (n > 127)
+            {
+              n = 127;
+            }
+        }
+    }
+  else if (n < -126)
+    {
+      y *= FLT_MIN * 0x1p24f;
+      n += 126 - 24;
+      if (n < -126)
+        {
+          y *= FLT_MIN * 0x1p24f;
+          n += 126 - 24;
+          if (n < -126)
+            {
+              n = -126;
+            }
+        }
+    }
+
+  u.i = (uint32_t)(0x7f + n) << 23;
+  x = y * u.f;
+  return x;
+}
+

libs/libc/math/lib_scalbnl.c

@@ -0,0 +1,144 @@
+/****************************************************************************
+ * libs/libc/math/lib_scalbnl.c
+ * get a long double number of x*2^n
+ *
+ * This file is copy from musl libc
+ * musl is an implementation of the C standard library built on top of the
+ * Linux system call API, including interfaces defined in the base language
+ * standard, POSIX, and widely agreed-upon extensions.
+ * musl is lightweight, fast, simple, free, and strives to be correct in
+ * the sense of standards-conformance and safety.
+ *
+ *
+ ****************************************************************************/
+
+/****************************************************************************
+ * Included Files
+ ****************************************************************************/
+
+#include <math.h>
+#include <stdint.h>
+#include <endian.h>
+#include <float.h>
+
+/****************************************************************************
+ * Pre-processor definitions
+ ****************************************************************************/
+
+#ifdef CONFIG_HAVE_LONG_DOUBLE
+
+#define LDOUBLE_MIN (0x1p-16382l)
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN
+union ldshape
+{
+  long double f;
+  struct
+    {
+      uint64_t m;
+      uint16_t se;
+    }i;
+};
+#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __BIG_ENDIAN
+
+/* This is the m68k variant of 80-bit long double,
+ * and this definition only works on archs where
+ * the alignment requirement of uint64_t is <= 4.
+ */
+
+union ldshape
+  {
+    long double f;
+    struct
+      {
+        uint16_t se;
+        uint16_t pad;
+        uint64_t m;
+      }i;
+};
+#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN
+union ldshape
+  {
+    long double f;
+    struct
+      {
+        uint64_t lo;
+        uint32_t mid;
+        uint16_t top;
+        uint16_t se;
+      }i;
+};
+#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __BIG_ENDIAN
+union ldshape
+  {
+    long double f;
+    struct
+      {
+        uint16_t se;
+        uint16_t top;
+        uint32_t mid;
+        uint64_t lo;
+      }i;
+};
+#else
+#error Unsupported long double representation
+#endif
+
+/****************************************************************************
+ * Public Functions
+ ****************************************************************************/
+
+/****************************************************************************
+ * Name: scalbnl
+ *
+ * Description:
+ *   get a long double number of x*2^n
+ *
+ ****************************************************************************/
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double scalbnl(long double x, int n)
+{
+  return scalbn(x, n);
+}
+#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
+long double scalbnl(long double x, int n)
+{
+  union ldshape u;
+
+  if (n > 16383)
+    {
+      x *= 0x1p16383l;
+      n -= 16383;
+      if (n > 16383)
+        {
+          x *= 0x1p16383l;
+          n -= 16383;
+          if (n > 16383)
+            {
+              n = 16383;
+            }
+        }
+    }
+  else if (n < -16382)
+    {
+      x *= LDOUBLE_MIN * 0x1p113l;
+      n += 16382 - 113;
+      if (n < -16382)
+        {
+          x *= LDOUBLE_MIN * 0x1p113l;
+          n += 16382 - 113;
+          if (n < -16382)
+            {
+              n = -16382;
+            }
+        }
+    }
+
+  u.f = 1.0;
+  u.i.se = 0x3fff + n;
+  return x * u.f;
+}
+#endif  /* LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 */
+#endif  /* CONFIG_HAVE_LONG_DOUBLE */
+