diff --git a/libs/libc/math/Make.defs b/libs/libc/math/Make.defs index 4b2cdfe7019..aa89ecdb2e5 100644 --- a/libs/libc/math/Make.defs +++ b/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 diff --git a/libs/libc/math/lib_scalbn.c b/libs/libc/math/lib_scalbn.c new file mode 100755 index 00000000000..93e2a060b99 --- /dev/null +++ b/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 +#include + +/**************************************************************************** + * 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 */ diff --git a/libs/libc/math/lib_scalbnf.c b/libs/libc/math/lib_scalbnf.c new file mode 100755 index 00000000000..7710a73c34d --- /dev/null +++ b/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 +#include + +/**************************************************************************** + * 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; +} + diff --git a/libs/libc/math/lib_scalbnl.c b/libs/libc/math/lib_scalbnl.c new file mode 100644 index 00000000000..0d93264537a --- /dev/null +++ b/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 +#include +#include +#include + +/**************************************************************************** + * 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 */ +