From 849d8d8d53cd514ffa63094e2e3c89327b7f82de Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miko=C5=82aj=20Mata=C5=82owski?= Date: Fri, 17 Jul 2026 17:25:38 +0200 Subject: [PATCH 1/3] libm: Move math library implementation, split implementation to seperate files YT: RTOS-1132 --- Makefile | 7 +- include/arch/aarch64/arch.h | 40 -- include/arch/armv7a/arch.h | 15 +- include/arch/armv8r/arch.h | 11 - include/arch/ia32/arch.h | 28 -- include/complex.h | 255 ------------- include/math.h | 230 ------------ include/math/consts.h | 71 ---- libm/Makefile | 68 ++++ libm/libmcs | 1 + {math => libm/phoenix}/common.c | 0 {math => libm/phoenix}/common.h | 0 libm/phoenix/compatibility.c | 49 +++ {math => libm/phoenix}/complex.c | 2 +- libm/phoenix/mathd/acosd.c | 51 +++ libm/phoenix/mathd/asind.c | 31 ++ libm/phoenix/mathd/atan2d.c | 73 ++++ libm/phoenix/mathd/atand.c | 58 +++ libm/phoenix/mathd/ceild.c | 35 ++ libm/phoenix/mathd/cosd.c | 67 ++++ libm/phoenix/mathd/coshd.c | 39 ++ libm/phoenix/mathd/expd.c | 58 +++ libm/phoenix/mathd/fabsd.c | 30 ++ libm/phoenix/mathd/floord.c | 34 ++ libm/phoenix/mathd/fmodd.c | 42 +++ libm/phoenix/mathd/frexpd.c | 45 +++ libm/phoenix/mathd/ldexpd.c | 55 +++ libm/phoenix/mathd/log10d.c | 23 ++ libm/phoenix/mathd/log2d.c | 22 ++ libm/phoenix/mathd/logd.c | 63 ++++ libm/phoenix/mathd/modfd.c | 57 +++ math/power.c => libm/phoenix/mathd/powd.c | 86 +---- libm/phoenix/mathd/roundd.c | 37 ++ libm/phoenix/mathd/sind.c | 70 ++++ math/hyper.c => libm/phoenix/mathd/sinhd.c | 47 +-- libm/phoenix/mathd/sqrtd.c | 69 ++++ libm/phoenix/mathd/tand.c | 41 ++ libm/phoenix/mathd/tanhd.c | 34 ++ libm/phoenix/mathd/truncd.c | 30 ++ libm/phoenix/mathf/acosf.c | 22 ++ libm/phoenix/mathf/asinf.c | 22 ++ libm/phoenix/mathf/atan2f.c | 22 ++ libm/phoenix/mathf/atanf.c | 22 ++ libm/phoenix/mathf/ceilf.c | 22 ++ libm/phoenix/mathf/cosf.c | 22 ++ libm/phoenix/mathf/coshf.c | 22 ++ libm/phoenix/mathf/expf.c | 22 ++ libm/phoenix/mathf/fabsf.c | 22 ++ libm/phoenix/mathf/floorf.c | 22 ++ libm/phoenix/mathf/fmodf.c | 22 ++ libm/phoenix/mathf/frexpf.c | 22 ++ libm/phoenix/mathf/ldexpf.c | 22 ++ libm/phoenix/mathf/log10f.c | 22 ++ libm/phoenix/mathf/logf.c | 22 ++ libm/phoenix/mathf/modff.c | 27 ++ libm/phoenix/mathf/powf.c | 22 ++ libm/phoenix/mathf/roundf.c | 22 ++ libm/phoenix/mathf/sinf.c | 21 ++ libm/phoenix/mathf/sinhf.c | 22 ++ libm/phoenix/mathf/sqrtf.c | 22 ++ libm/phoenix/mathf/tanf.c | 22 ++ libm/phoenix/mathf/tanhf.c | 22 ++ libm/phoenix/mathf/truncf.c | 22 ++ math/Makefile | 7 - math/exp.c | 411 --------------------- math/trig.c | 293 --------------- 66 files changed, 1727 insertions(+), 1490 deletions(-) delete mode 100644 include/complex.h delete mode 100644 include/math.h delete mode 100644 include/math/consts.h create mode 100644 libm/Makefile create mode 160000 libm/libmcs rename {math => libm/phoenix}/common.c (100%) rename {math => libm/phoenix}/common.h (100%) create mode 100644 libm/phoenix/compatibility.c rename {math => libm/phoenix}/complex.c (96%) create mode 100644 libm/phoenix/mathd/acosd.c create mode 100644 libm/phoenix/mathd/asind.c create mode 100644 libm/phoenix/mathd/atan2d.c create mode 100644 libm/phoenix/mathd/atand.c create mode 100644 libm/phoenix/mathd/ceild.c create mode 100644 libm/phoenix/mathd/cosd.c create mode 100644 libm/phoenix/mathd/coshd.c create mode 100644 libm/phoenix/mathd/expd.c create mode 100644 libm/phoenix/mathd/fabsd.c create mode 100644 libm/phoenix/mathd/floord.c create mode 100644 libm/phoenix/mathd/fmodd.c create mode 100644 libm/phoenix/mathd/frexpd.c create mode 100644 libm/phoenix/mathd/ldexpd.c create mode 100644 libm/phoenix/mathd/log10d.c create mode 100644 libm/phoenix/mathd/log2d.c create mode 100644 libm/phoenix/mathd/logd.c create mode 100644 libm/phoenix/mathd/modfd.c rename math/power.c => libm/phoenix/mathd/powd.c (56%) create mode 100644 libm/phoenix/mathd/roundd.c create mode 100644 libm/phoenix/mathd/sind.c rename math/hyper.c => libm/phoenix/mathd/sinhd.c (60%) create mode 100644 libm/phoenix/mathd/sqrtd.c create mode 100644 libm/phoenix/mathd/tand.c create mode 100644 libm/phoenix/mathd/tanhd.c create mode 100644 libm/phoenix/mathd/truncd.c create mode 100644 libm/phoenix/mathf/acosf.c create mode 100644 libm/phoenix/mathf/asinf.c create mode 100644 libm/phoenix/mathf/atan2f.c create mode 100644 libm/phoenix/mathf/atanf.c create mode 100644 libm/phoenix/mathf/ceilf.c create mode 100644 libm/phoenix/mathf/cosf.c create mode 100644 libm/phoenix/mathf/coshf.c create mode 100644 libm/phoenix/mathf/expf.c create mode 100644 libm/phoenix/mathf/fabsf.c create mode 100644 libm/phoenix/mathf/floorf.c create mode 100644 libm/phoenix/mathf/fmodf.c create mode 100644 libm/phoenix/mathf/frexpf.c create mode 100644 libm/phoenix/mathf/ldexpf.c create mode 100644 libm/phoenix/mathf/log10f.c create mode 100644 libm/phoenix/mathf/logf.c create mode 100644 libm/phoenix/mathf/modff.c create mode 100644 libm/phoenix/mathf/powf.c create mode 100644 libm/phoenix/mathf/roundf.c create mode 100644 libm/phoenix/mathf/sinf.c create mode 100644 libm/phoenix/mathf/sinhf.c create mode 100644 libm/phoenix/mathf/sqrtf.c create mode 100644 libm/phoenix/mathf/tanf.c create mode 100644 libm/phoenix/mathf/tanhf.c create mode 100644 libm/phoenix/mathf/truncf.c delete mode 100644 math/Makefile delete mode 100644 math/exp.c delete mode 100644 math/trig.c diff --git a/Makefile b/Makefile index 126227e0..506c346a 100644 --- a/Makefile +++ b/Makefile @@ -50,7 +50,7 @@ include ctype/Makefile include err/Makefile include errno/Makefile include locale/Makefile -include math/Makefile +include libm/Makefile include misc/Makefile include net/Makefile include netinet/Makefile @@ -86,7 +86,8 @@ install: install-headers install-libs install-headers: $(SRCHEADERS) @echo INSTALL "$(HEADERS_INSTALL_DIR)/*"; \ mkdir -p "$(HEADERS_INSTALL_DIR)"; \ - cp -a include/* "$(HEADERS_INSTALL_DIR)"; + cp -a include/* "$(HEADERS_INSTALL_DIR)"; \ + cp -a libm/libmcs/libm/include/* "$(HEADERS_INSTALL_DIR)"; # TODO: remove `rm crt0.o` when we will be sure it's not a symlink to libphoenix.a anymore install-libs: $(LIB_TARGETS) @@ -108,4 +109,4 @@ clean: ifneq ($(filter clean,$(MAKECMDGOALS)),) $(shell rm -rf $(BUILD_DIR)) $(shell rm -f string/*.inc) -endif +endif \ No newline at end of file diff --git a/include/arch/aarch64/arch.h b/include/arch/aarch64/arch.h index 5ed86bc3..5aaf38c5 100644 --- a/include/arch/aarch64/arch.h +++ b/include/arch/aarch64/arch.h @@ -23,46 +23,6 @@ #define __MEMSET #define __MEMMOVE -#ifndef __SOFTFP__ -/* clang-format off */ -#define __IEEE754_SQRT -#define __ieee754_sqrt(x) ({ double a = (x); __asm__ ("fsqrt %d0, %d1" : "=w"(a) : "w"(a)); a; }) - -#define __IEEE754_SQRTF -#define __ieee754_sqrtf(x) ({ float a = (x); __asm__ ("fsqrt %s0, %s1" : "=w"(a) : "w"(a)); a; }) - -#define __IEEE754_FABS -#define __ieee754_fabs(x) ({ double a = (x); __asm__ ("fabs %d0, %d1" : "=w"(a) : "w"(a)); a; }) - -#define __IEEE754_FABSF -#define __ieee754_fabsf(x) ({ float a = (x); __asm__ ("fabs %s0, %s1" : "=w"(a) : "w"(a)); a; }) - -#define __IEEE754_CEIL -#define __ieee754_ceil(x) ({ double a = (x); __asm__ ("frintp %d0, %d1" : "=w"(a) : "w"(a)); a; }) - -#define __IEEE754_CEILF -#define __ieee754_ceilf(x) ({ float a = (x); __asm__ ("frintp %s0, %s1" : "=w"(a) : "w"(a)); a; }) - -#define __IEEE754_FLOOR -#define __ieee754_floor(x) ({ double a = (x); __asm__ ("frintm %d0, %d1" : "=w"(a) : "w"(a)); a; }) - -#define __IEEE754_FLOORF -#define __ieee754_floorf(x) ({ float a = (x); __asm__ ("frintm %s0, %s1" : "=w"(a) : "w"(a)); a; }) - -#define __IEEE754_ROUND -#define __ieee754_round(x) ({ double a = (x); __asm__ ("frinta %d0, %d1" : "=w"(a) : "w"(a)); a; }) - -#define __IEEE754_ROUNDF -#define __ieee754_roundf(x) ({ float a = (x); __asm__ ("frinta %s0, %s1" : "=w"(a) : "w"(a)); a; }) - -#define __IEEE754_TRUNC -#define __ieee754_trunc(x) ({ double a = (x); __asm__ ("frintz %d0, %d1" : "=w"(a) : "w"(a)); a; }) - -#define __IEEE754_TRUNCF -#define __ieee754_truncf(x) ({ float a = (x); __asm__ ("frintz %s0, %s1" : "=w"(a) : "w"(a)); a; }) -/* clang-format on */ -#endif - #define __LIBPHOENIX_ARCH_TLS_SUPPORTED #endif diff --git a/include/arch/armv7a/arch.h b/include/arch/armv7a/arch.h index cb2ed681..a6982ff5 100644 --- a/include/arch/armv7a/arch.h +++ b/include/arch/armv7a/arch.h @@ -16,8 +16,8 @@ #ifndef _LIBPHOENIX_ARCH_ARMV7A_ARCH_H_ #define _LIBPHOENIX_ARCH_ARMV7A_ARCH_H_ -#define __ARCH_STDINT -#define __ARCH_LIMITS +#define __ARCH_STDINT +#define __ARCH_LIMITS #define __MEMCPY #define __MEMCMP @@ -30,17 +30,6 @@ #define __STRNCPY #define __MEMMOVE - -#if defined(__ARM_PCS_VFP) || (defined(__VFP_FP__) && !defined(__SOFTFP__)) -#if defined(__ARM_FP) && (__ARM_FP & 8) != 0 -#define __IEEE754_SQRT -#define __ieee754_sqrt(x) ({ double a = (x); __asm__ volatile ("vsqrt.f64 %P0, %P1" : "=w"(a) : "w"(a)); a; }) -#endif - -#define __IEEE754_SQRTF -#define __ieee754_sqrtf(x) ({ float a = (x); __asm__ volatile ("vsqrt.f32 %0, %1" : "=t"(a) : "t"(a)); a; }) -#endif - #define __LIBPHOENIX_ARCH_TLS_SUPPORTED #endif diff --git a/include/arch/armv8r/arch.h b/include/arch/armv8r/arch.h index 2831c1f9..2c4aad22 100644 --- a/include/arch/armv8r/arch.h +++ b/include/arch/armv8r/arch.h @@ -30,17 +30,6 @@ #define __STRNCPY #define __MEMMOVE - -#if defined(__ARM_PCS_VFP) || (defined(__VFP_FP__) && !defined(__SOFTFP__)) -#if defined(__ARM_FP) && (__ARM_FP & 8) != 0 -#define __IEEE754_SQRT -#define __ieee754_sqrt(x) ({ double a = (x); __asm__ volatile ("vsqrt.f64 %P0, %P1" : "=w"(a) : "w"(a)); a; }) -#endif - -#define __IEEE754_SQRTF -#define __ieee754_sqrtf(x) ({ float a = (x); __asm__ volatile ("vsqrt.f32 %0, %1" : "=t"(a) : "t"(a)); a; }) -#endif - #define __LIBPHOENIX_ARCH_TLS_SUPPORTED #endif diff --git a/include/arch/ia32/arch.h b/include/arch/ia32/arch.h index c00dad97..b05a309d 100644 --- a/include/arch/ia32/arch.h +++ b/include/arch/ia32/arch.h @@ -23,34 +23,6 @@ #define __MEMCPY #define __MEMSET -#ifndef __SOFTFP__ -#define __IEEE754_SQRT - -static inline double __ieee754_sqrt(double x) -{ - unsigned short newcw, savecw; - unsigned short register reg; - - /* clang-format off */ - __asm__ volatile ( - "fstcw %w[savecw]\n\t" /* save control word */ - "mov %w[savecw], %w[reg]\n\t" /* reg <- (memptr) */ - "and $0xf0ff, %w[reg]\n\t" /* inherit exception mask and infinity */ - "or $0x200, %w[reg]\n\t" /* set double precision, round to nearest-even */ - "mov %w[reg], %w[newcw]\n\t" /* (memptr) <- reg */ - "fldcw %w[newcw]\n\t" /* and apply (memptr) new control word */ - "fsqrt\n\t" /* calculate double precision sqrt */ - "fldcw %w[savecw]" /* restore control word */ - : "+t" (x), [reg] "=&r" (reg) - : [newcw] "m" (newcw), [savecw] "m" (savecw) - : "cc", "memory"); - /* clang-format on */ - - return x; -} - -#endif - #define __LIBPHOENIX_ARCH_TLS_SUPPORTED #endif diff --git a/include/complex.h b/include/complex.h deleted file mode 100644 index 9c7a3427..00000000 --- a/include/complex.h +++ /dev/null @@ -1,255 +0,0 @@ -/* - * Phoenix-RTOS - * - * libphoenix - * - * mathematical functions on complex numbers - * - * Copyright 2023 Phoenix Systems - * Author: Gerard Swiderski - * - * This file is part of Phoenix-RTOS. - * - * %LICENSE% - */ - -#ifndef _COMPLEX_H_ -#define _COMPLEX_H_ - -/* define complex numbers specifers */ -#define complex _Complex -#define _Complex_I 1.0fi -#define I _Complex_I - - -#ifdef __cplusplus -extern "C" { -#endif - - -static inline double creal(double complex z) -{ - return (__real__(z)); -} - - -static inline float crealf(float complex z) -{ - return (__real__(z)); -} - - -static inline double cimag(double complex z) -{ - return (__imag__(z)); -} - - -static inline float cimagf(float complex z) -{ - return (__imag__(z)); -} - - -static inline double complex conj(double complex z) -{ - return __builtin_complex((__real__(z)), -(__imag__(z))); -} - - -static inline float complex conjf(float complex z) -{ - return __builtin_complex((__real__(z)), -(__imag__(z))); -} - - -double cabs(double complex); - - -float cabsf(float complex); - - -double carg(double complex); - - -float cargf(float complex); - - -double complex cexp(double complex); - - -float complex cexpf(float complex); - - -/* - * TODO: functions not yet implemented - */ - -long double cabsl(long double complex); - - -double complex cacos(double complex); - - -float complex cacosf(float complex); - - -double complex cacosh(double complex); - - -float complex cacoshf(float complex); - - -long double complex cacoshl(long double complex); - - -long double complex cacosl(long double complex); - - -long double cargl(long double complex); - - -double complex casin(double complex); - - -float complex casinf(float complex); - - -double complex casinh(double complex); - - -float complex casinhf(float complex); - - -long double complex casinhl(long double complex); - - -long double complex casinl(long double complex); - - -double complex catan(double complex); - - -float complex catanf(float complex); - - -double complex catanh(double complex); - - -float complex catanhf(float complex); - - -long double complex catanhl(long double complex); - - -long double complex catanl(long double complex); - - -double complex ccos(double complex); - - -float complex ccosf(float complex); - - -double complex ccosh(double complex); - - -float complex ccoshf(float complex); - - -long double complex ccoshl(long double complex); - - -long double complex ccosl(long double complex); - - -long double complex cexpl(long double complex); - - -long double cimagl(long double complex); - - -double complex clog(double complex); - - -float complex clogf(float complex); - - -long double complex clogl(long double complex); - - -long double complex conjl(long double complex); - - -double complex cpow(double complex, double complex); - - -float complex cpowf(float complex, float complex); - - -long double complex cpowl(long double complex, long double complex); - - -double complex cproj(double complex); - - -float complex cprojf(float complex); - - -long double complex cprojl(long double complex); - - -long double creall(long double complex); - - -double complex csin(double complex); - - -float complex csinf(float complex); - - -double complex csinh(double complex); - - -float complex csinhf(float complex); - - -long double complex csinhl(long double complex); - - -long double complex csinl(long double complex); - - -double complex csqrt(double complex); - - -float complex csqrtf(float complex); - - -long double complex csqrtl(long double complex); - - -double complex ctan(double complex); - - -float complex ctanf(float complex); - - -double complex ctanh(double complex); - - -float complex ctanhf(float complex); - - -long double complex ctanhl(long double complex); - - -long double complex ctanl(long double complex); - - -#ifdef __cplusplus -} -#endif - - -#endif /* end of _COMPLEX_H_ */ diff --git a/include/math.h b/include/math.h deleted file mode 100644 index 3353192c..00000000 --- a/include/math.h +++ /dev/null @@ -1,230 +0,0 @@ -/* - * Phoenix-RTOS - * - * libphoenix - * - * math.h - * - * Copyright 2017, 2022 Phoenix Systems - * Author: Aleksander Kaminski, Damian Loewnau - * - * This file is part of Phoenix-RTOS. - * - * %LICENSE% - */ - -#ifndef _LIBPHOENIX_MATH_H_ -#define _LIBPHOENIX_MATH_H_ - - -#include -#include - - -#ifdef __cplusplus -extern "C" { -#endif - -#define MATH_ERRNO 1 -#define MATH_ERREXCEPT 2 -#define math_errhandling (MATH_ERRNO) - -#define FP_NAN 0 -#define FP_INFINITE 1 -#define FP_ZERO 2 -#define FP_SUBNORMAL 3 -#define FP_NORMAL 4 - - -#define fpclassify(x) __builtin_fpclassify(FP_NAN, FP_INFINITE, FP_NORMAL, FP_SUBNORMAL, FP_ZERO, x) -#define isfinite(x) __builtin_isfinite(x) -#define isgreater(x, y) __builtin_isgreater(x, y) -#define isgreaterequal(x, y) __builtin_isgreaterequal(x, y) -#define isinf(x) __builtin_isinf(x) -#define isless(x, y) __builtin_isless(x, y) -#define islessequal(x, y) __builtin_islessequal(x, y) -#define islessgreater(x, y) __builtin_islessgreater(x, y) -#define isnan(x) __builtin_isnan(x) -#define signbit(x) __builtin_signbit(x) -#define isnormal(x) __builtin_isnormal(x) -#define isunordered(x, y) __builtin_isunordered(x, y) - -#define HUGE_VAL __builtin_huge_val() -#define HUGE_VALF __builtin_huge_valf() -#define HUGE_VALL __builtin_huge_vall() -#define INFINITY __builtin_inff() -#define NAN __builtin_nanf("") - - -/* Trigonometric functions */ - - -/* Returns the cosine of an angle of x radians. */ -extern double cos(double x); - - -/* Returns the sine of an angle of x radians. */ -extern double sin(double x); - - -/* Returns the tangent of an angle of x radians. */ -extern double tan(double x); - - -/* Returns the principal value of the arc cosine of x, expressed in radians. */ -extern double acos(double x); - - -/* Returns the principal value of the arc sine of x, expressed in radians. */ -extern double asin(double x); - - -/* Returns the principal value of the arc tan of x, expressed in radians. */ -extern double atan(double x); - - -/* Returns the principal value of the arc tangent of y/x, expressed in radians. */ -extern double atan2(double y, double x); - - -/* Hyperbolic functions */ - - -/* Returns the hyperbolic cosine of x. */ -extern double cosh(double x); - - -/* Returns the hyperbolic sine of x. */ -extern double sinh(double x); - - -/* Returns the hyperbolic tangent of x. */ -extern double tanh(double x); - - -/* Exponential and logarithmic functions */ - - -/* Returns the base-e exponential function of x, which is e raised to the power x: e^x. */ -extern double exp(double x); - - -/* Breaks the floating point number x into its binary significand - * (a floating point with an absolute value between 0.5(included) and 1.0(excluded)) - * and an integral exponent for 2. */ -extern double frexp(double x, int *exp); - - -/* Returns the result of multiplying x (the significand) by 2 raised to the power of exp (the exponent). */ -extern double ldexp(double x, int exp); - - -/* Returns the natural logarithm of x. */ -extern double log(double x); - - -/* Returns the common (base-2) logarithm of x. */ -extern double log2(double x); - - -/* Returns the common (base-10) logarithm of x. */ -extern double log10(double x); - - -/* Breaks x into an integral and a fractional part. */ -extern double modf(double x, double *intpart); -extern float modff(float x, float *intpart); - - -/* Power functions */ - - -/* Returns base raised to the power exponent. */ -extern double pow(double base, double exponent); - - -/* Returns the square root of x. */ -extern double sqrt(double x); - - -/* Rounding and remainder functions */ - - -/* Rounds x upward, returning the smallest integral value that is not less than x. */ -extern double ceil(double x); -extern float ceilf(float x); - - -/* Rounds x downward, returning the largest integral value that is not greater than x. */ -extern double floor(double x); -extern float floorf(float x); - - -/* Returns the floating-point remainder of numer/denom (rounded towards zero). */ -extern double fmod(double numer, double denom); - - -/* Return the integral value nearest to x */ -extern double round(double x); -extern float roundf(float x); - - -/* Rounds x toward zero, returning the nearest integral value that is not larger in magnitude than x. */ -extern double trunc(double x); -extern float truncf(float x); - - -/* Miscellaneous */ - - -/* Returns the absolute value of x: |x|. */ -extern double fabs(double x); -extern float fabsf(float x); - - -/* C99 extensions */ -float cosf(float x); -float sinf(float x); -float tanf(float x); -float acosf(float x); -float asinf(float x); -float atanf(float x); -float atan2f(float y, float x); -float coshf(float x); -float sinhf(float x); -float tanhf(float x); -float expf(float x); -float frexpf(float x, int *exp); -float ldexpf(float x, int exp); -float logf(float x); -float log10f(float x); -float powf(float base, float exponent); -float sqrtf(float x); -float fmodf(float num, float denom); - - -#define cosf(x) ((float)cos(x)) -#define sinf(x) ((float)sin(x)) -#define tanf(x) ((float)tan(x)) -#define acosf(x) ((float)acos(x)) -#define asinf(x) ((float)asin(x)) -#define atanf(x) ((float)atan(x)) -#define atan2f(y, x) ((float)atan2(y, x)) -#define coshf(x) ((float)cosh(x)) -#define sinhf(x) ((float)sinh(x)) -#define tanhf(x) ((float)tanh(x)) -#define expf(x) ((float)exp(x)) -#define frexpf(x, exp) ((float)frexp(x, exp)) -#define ldexpf(x, exp) ((float)ldexp(x, exp)) -#define logf(x) ((float)log(x)) -#define log10f(x) ((float)log10(x)) -#define powf(base, exponent) ((float)pow(base, exponent)) -#define fmodf(num, denom) ((float)fmod(num, denom)) - - -#ifdef __cplusplus -} -#endif - - -#endif diff --git a/include/math/consts.h b/include/math/consts.h deleted file mode 100644 index 40ac3e7d..00000000 --- a/include/math/consts.h +++ /dev/null @@ -1,71 +0,0 @@ -/* - * Phoenix-RTOS - * - * libphoenix - * - * math.h constants - * - * Copyright 2017 Phoenix Systems - * Author: Aleksander Kaminski - * - * This file is part of Phoenix-RTOS. - * - * %LICENSE% - */ - -#ifndef _LIBPHOENIX_MATH_CONSTS_H_ -#define _LIBPHOENIX_MATH_CONSTS_H_ - -/* The base of natural logarithms. */ -#define M_E 2.7182818284590452354 - - -/* The logarithm to base 2 of M_E. */ -#define M_LOG2E 1.4426950408889634074 - - -/* The logarithm to base 10 of M_E. */ -#define M_LOG10E 0.43429448190325182765 - - -/* The natural logarithm of 2. */ -#define M_LN2 0.69314718055994530942 - - -/* The natural logarithm of 10. */ -#define M_LN10 2.30258509299404568402 - - -/* Pi, the ratio of a circle's circumference to its diameter. */ -#define M_PI 3.14159265358979323846 - - -/* Pi divided by two. */ -#define M_PI_2 1.57079632679489661923 - - -/* Pi divided by four. */ -#define M_PI_4 0.78539816339744830962 - - -/* The reciprocal of pi (1/pi) */ -#define M_1_PI 0.31830988618379067154 - - -/* Two times the reciprocal of pi. */ -#define M_2_PI 0.63661977236758134308 - - -/* Two times the reciprocal of the square root of pi. * */ -#define M_2_SQRTPI 1.12837916709551257390 - - -/* The square root of two. */ -#define M_SQRT2 1.41421356237309504880 - - -/* The reciprocal of the square root of two (also the square root of 1/2). */ -#define M_SQRT1_2 0.70710678118654752440 - - -#endif diff --git a/libm/Makefile b/libm/Makefile new file mode 100644 index 00000000..40a8dad1 --- /dev/null +++ b/libm/Makefile @@ -0,0 +1,68 @@ +# +# Makefile for math library +# +# Copyright 2025 Phoenix Systems +# Author: Mikolaj Matalowski +# +# This file is part of Phoenix-RTOS. +# +# %LICENSE% +# + +LIBM_USE_LIBMCS ?= n +LIBM_LIBMCS_DAZ ?= n +LIBM_WANT_COMPLEX ?= y +LIBM_USE_HW ?= y + +# Set options for denormals on the FPU +ifeq ($(LIBM_LIBMCS_DAZ), y) + CPPFLAGS += -DLIBMCS_FPU_DAZ +endif + +ifeq ($(LIBM_WANT_COMPLEX), y) + CPPFLAGS += -DLIBMCS_WANT_COMPLEX +endif + +CFLAGS += -Ilibm/libmcs/libm/include + +# Add source files for hardware accelerated math +HW_OBJS := +ifeq ($(LIBM_USE_HW), y) + include libm/arch/Makefile +endif + +SW_OBJS := +# Build for libmcs +ifeq ($(LIBM_USE_LIBMCS), y) + ifeq ($(LIBM_WANT_COMPLEX), y) + CPPFLAGS += -DLIBM_WANT_COMPLEX + SW_OBJS += $(addprefix $(PREFIX_O), $(patsubst %.c,%.o,$(wildcard libm/libmcs/libm/complexd/internal/*.c))) + SW_OBJS += $(addprefix $(PREFIX_O), $(patsubst %.c,%.o,$(wildcard libm/libmcs/libm/complexd/*.c))) + SW_OBJS += $(addprefix $(PREFIX_O), $(patsubst %.c,%.o,$(wildcard libm/libmcs/libm/complexf/internal/*.c))) + SW_OBJS += $(addprefix $(PREFIX_O), $(patsubst %.c,%.o,$(wildcard libm/libmcs/libm/complexf/*.c))) + endif + + # fenv.c implementation is target specific + SW_OBJS += $(addprefix $(PREFIX_O), $(patsubst %.c,%.o, $(filter-out %fenv.c,$(wildcard libm/libmcs/libm/common/*.c)))) + SW_OBJS += $(addprefix $(PREFIX_O), $(patsubst %.c,%.o,$(wildcard libm/libmcs/libm/mathf/*.c))) + SW_OBJS += $(addprefix $(PREFIX_O), $(patsubst %.c,%.o,$(wildcard libm/libmcs/libm/mathf/internal/*.c))) + SW_OBJS += $(addprefix $(PREFIX_O), $(patsubst %.c,%.o,$(wildcard libm/libmcs/libm/mathd/*.c))) + SW_OBJS += $(addprefix $(PREFIX_O), $(patsubst %.c,%.o,$(wildcard libm/libmcs/libm/mathd/internal/*.c))) +else + ifeq ($(LIBM_WANT_COMPLEX), y) + SW_OBJS += $(addprefix $(PREFIX_O), $(patsubst %.c,%.o, libm/phoenix/complex.c)) + endif + + # Build for phoenix implementation + SW_OBJS += $(addprefix $(PREFIX_O), $(patsubst %.c,%.o, $(patsubst %complex.c,, $(wildcard libm/phoenix/*.c)))) + SW_OBJS += $(addprefix $(PREFIX_O), $(patsubst %.c,%.o,$(wildcard libm/phoenix/mathd/*.c))) + SW_OBJS += $(addprefix $(PREFIX_O), $(patsubst %.c,%.o,$(wildcard libm/phoenix/mathf/*.c))) + + CFLAGS += -Ilibm/phoenix + # This part is needed for libmcs headers compatibility + CFLAGS += -include libm/libmcs/libm/include/internal_config.h +endif + +# Clear functions that already implemented in hardware +OBJS += $(foreach f,$(SW_OBJS),$(if $(filter $(notdir $(f)),$(notdir $(HW_OBJS))),,$(f))) +OBJS += $(HW_OBJS) \ No newline at end of file diff --git a/libm/libmcs b/libm/libmcs new file mode 160000 index 00000000..598059c4 --- /dev/null +++ b/libm/libmcs @@ -0,0 +1 @@ +Subproject commit 598059c4ad6e76ad90454c010de86f08a1d0d49d diff --git a/math/common.c b/libm/phoenix/common.c similarity index 100% rename from math/common.c rename to libm/phoenix/common.c diff --git a/math/common.h b/libm/phoenix/common.h similarity index 100% rename from math/common.h rename to libm/phoenix/common.h diff --git a/libm/phoenix/compatibility.c b/libm/phoenix/compatibility.c new file mode 100644 index 00000000..6350d9fe --- /dev/null +++ b/libm/phoenix/compatibility.c @@ -0,0 +1,49 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * Compatibility source file for libmcs header file + * + * Copyright 2025 Phoenix Systems + * Author: Mikolaj Matalowski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ +#include + + +int __fpclassifyf(float x) +{ + return __builtin_fpclassify(FP_NAN, FP_INFINITE, FP_NORMAL, FP_SUBNORMAL, FP_ZERO, x); +} + + +int __fpclassifyd(double x) +{ + return __builtin_fpclassify(FP_NAN, FP_INFINITE, FP_NORMAL, FP_SUBNORMAL, FP_ZERO, x); +} + + +float nanf(const char *) +{ + return __builtin_nanf(""); +} + + +int __signbitf(float x) +{ + return __builtin_signbitf(x); +} + + +int __signbitd(double x) +{ + return __builtin_signbit(x); +} + + +const float __inff = __builtin_inff(); +const double __infd = __builtin_inf(); diff --git a/math/complex.c b/libm/phoenix/complex.c similarity index 96% rename from math/complex.c rename to libm/phoenix/complex.c index c12d4e60..2976b897 100644 --- a/math/complex.c +++ b/libm/phoenix/complex.c @@ -10,7 +10,7 @@ * * This file is part of Phoenix-RTOS. * - * %LICENSE% + * SPDX-License-Identifier: BSD-3-Clause */ #include diff --git a/libm/phoenix/mathd/acosd.c b/libm/phoenix/mathd/acosd.c new file mode 100644 index 00000000..59532000 --- /dev/null +++ b/libm/phoenix/mathd/acosd.c @@ -0,0 +1,51 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * cos, sin, tan, acos, asin, atan + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Aleksander Kaminski, Jakub Smolaga + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include +#include + + +/* Calculates value of arc cosine using secant method */ +double acos(double x) +{ + double xa = 0, xb = M_PI, ya, yb, t; + int i; + + if (isnan(x) != 0) { + return NAN; + } + + if ((x > 1.0) || (x < -1.0)) { + errno = EDOM; + return NAN; + } + + for (i = 0; i < 16; ++i) { + ya = cos(xa) - x; + yb = cos(xb) - x; + + t = ya - yb; + + if (t == 0.0) { + break; + } + + t = ((ya * xb) - (yb * xa)) / t; + xa = xb; + xb = t; + } + + return xb; +} diff --git a/libm/phoenix/mathd/asind.c b/libm/phoenix/mathd/asind.c new file mode 100644 index 00000000..462deeca --- /dev/null +++ b/libm/phoenix/mathd/asind.c @@ -0,0 +1,31 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * cos, sin, tan, acos, asin, atan + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Aleksander Kaminski, Jakub Smolaga + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +/* Calculates value of arc sine using asin(x) = pi/2 - acos(x) relationship. */ +double asin(double x) +{ + if (isnan(x) != 0) { + return NAN; + } + + if (x == 0.0) { + return x; + } + + return M_PI_2 - acos(x); +} diff --git a/libm/phoenix/mathd/atan2d.c b/libm/phoenix/mathd/atan2d.c new file mode 100644 index 00000000..8aa7f354 --- /dev/null +++ b/libm/phoenix/mathd/atan2d.c @@ -0,0 +1,73 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * cos, sin, tan, acos, asin, atan + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Aleksander Kaminski, Jakub Smolaga + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +double atan2(double y, double x) +{ + if ((isnan(y) != 0) || (isnan(x) != 0)) { + return NAN; + } + + if (y == 0.0) { + if (x < 0.0 || signbit(x) != 0) { + // x is negative or -0.0 + return (signbit(y) != 0) ? -M_PI : M_PI; + } + else { + // x is positive or +0.0 + return y; + } + } + + if ((isinf(y) != 0) && (isinf(x) != 0)) { + int sy = signbit(y); + int sx = signbit(x); + + if ((sy == 0) && (sx != 0)) { + return M_PI_4 * 3; + } + else if ((sy != 0) && (sx != 0)) { + return -M_PI_4 * 3; + } + else if ((sy == 0) && (sx == 0)) { + return M_PI_4; + } + else if ((sy != 0) && (sx == 0)) { + return -M_PI_4; + } + } + + if (x > 0.0) { + return atan(y / x); + } + else if (x < 0.0) { + if (y >= 0.0) { + return atan(y / x) + M_PI; + } + + return atan(y / x) - M_PI; + } + + if (y > 0.0) { + return M_PI_2; + } + else if (y < 0.0) { + return -M_PI_2; + } + + return 0.0; +} diff --git a/libm/phoenix/mathd/atand.c b/libm/phoenix/mathd/atand.c new file mode 100644 index 00000000..73c3b42e --- /dev/null +++ b/libm/phoenix/mathd/atand.c @@ -0,0 +1,58 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * cos, sin, tan, acos, asin, atan + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Aleksander Kaminski, Jakub Smolaga + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +static const double atan_wi[] = { 0.2152638534631578, 0.2152638534631578, 0.2051984637212956, + 0.2051984637212956, 0.1855383974779378, 0.1855383974779378, 0.1572031671581935, + 0.1572031671581935, 0.1215185706879032, 0.1215185706879032, 0.0801580871597602, + 0.0801580871597602, 0.0351194603317519, 0.0351194603317519 }; + + +static const double atan_xi[] = { -0.1080549487073437, 0.1080549487073437, -0.3191123689278897, + 0.3191123689278897, -0.5152486363581541, 0.5152486363581541, -0.6872929048116855, + 0.6872929048116855, -0.8272013150697650, 0.8272013150697650, -0.9284348836635735, + 0.9284348836635735, -0.9862838086968123, 0.9862838086968123 }; + + +double atan(double x) +{ + double res = 1.0, h, a; + int i, s; + + if (isnan(x) != 0) { + return NAN; + } + + s = (x < 0.0) ? -1 : 1; + x *= (double)s; + h = x / 2; + + if (x == 0.0) { + return x; + } + + if (x > 1.0) { + return ((M_PI_2 - atan(1.0 / x)) * (double)s); + } + + for (i = 0, res = 0.0; i < sizeof(atan_wi) / sizeof(atan_wi[0]); ++i) { + a = (h * atan_xi[i]) + h; + res += atan_wi[i] / ((a * a) + 1); + } + + return (res * h * (double)s); +} diff --git a/libm/phoenix/mathd/ceild.c b/libm/phoenix/mathd/ceild.c new file mode 100644 index 00000000..0e56172f --- /dev/null +++ b/libm/phoenix/mathd/ceild.c @@ -0,0 +1,35 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + + +#include + + +double ceil(double x) +{ + double ipart, fpart; + + if (isnan(x) != 0) { + return NAN; + } + + fpart = modf(x, &ipart); + + if ((x > 0.0) && ((fpart + x) != x)) { + ipart += 1.0; + } + + return ipart; +} diff --git a/libm/phoenix/mathd/cosd.c b/libm/phoenix/mathd/cosd.c new file mode 100644 index 00000000..5d49b597 --- /dev/null +++ b/libm/phoenix/mathd/cosd.c @@ -0,0 +1,67 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * cos, sin, tan, acos, asin, atan + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Aleksander Kaminski, Jakub Smolaga + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include +#include + + +/* Calculates value of cosine using Maclaurin series via Horner's method. */ +double cos(double x) +{ + double res, xpow; + + if (isnan(x) != 0) { + return NAN; + } + + if (isinf(x) != 0) { + errno = EDOM; + return NAN; + } + + /* Normalize argument to -2*PI < x < 2*PI */ + x = fmod(x, 2.0 * M_PI); + + /* Normalize further to -PI < x < PI */ + if (x > M_PI) { + x -= 2.0 * M_PI; + } + else if (x < -M_PI) { + x += 2.0 * M_PI; + } + + if (x > M_PI_2) { + return -sin(x - M_PI_2); + } + else if (x < -M_PI_2) { + return sin(x + M_PI_2); + } + + xpow = x * x; + + res = 0x1.e542ba4020225p-62; /* 1/20! */ + res = res * xpow - 0x1.6827863b97d97p-53; /* 1/18! */ + res = res * xpow + 0x1.ae7f3e733b81fp-45; /* 1/16! */ + res = res * xpow - 0x1.93974a8c07c9dp-37; /* 1/14! */ + res = res * xpow + 0x1.1eed8eff8d898p-29; /* 1/12! */ + res = res * xpow - 0x1.27e4fb7789f5cp-22; /* 1/10! */ + res = res * xpow + 0x1.a01a01a01a01ap-16; /* 1/8! */ + res = res * xpow - 0x1.6c16c16c16c17p-10; /* 1/6! */ + res = res * xpow + 0x1.5555555555555p-5; /* 1/4! */ + res = res * xpow - 0x1p-1; /* 1/2! */ + res = res * xpow + 0x1p0; /* 1/0! */ + + return res; +} diff --git a/libm/phoenix/mathd/coshd.c b/libm/phoenix/mathd/coshd.c new file mode 100644 index 00000000..641a9357 --- /dev/null +++ b/libm/phoenix/mathd/coshd.c @@ -0,0 +1,39 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * hyberbolic functions + * + * Copyright 2018 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +double cosh(double x) +{ + double y; + + if (isnan(x) != 0) { + return NAN; + } + + if (isinf(x) != 0) { + return INFINITY; + } + + /* Make sure exp(x) is not infinity */ + if (x < 709.78) { + return ((exp(x) + exp(-x)) / 2.0); + } + else { + y = cosh(x / 2.0); + return ((2.0 * y * y) - 1.0); + } +} diff --git a/libm/phoenix/mathd/expd.c b/libm/phoenix/mathd/expd.c new file mode 100644 index 00000000..a79f9c7b --- /dev/null +++ b/libm/phoenix/mathd/expd.c @@ -0,0 +1,58 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include +#include +#include "common.h" + + +/* Uses quick powering and Maclaurin series to calculate value of e^x */ +double exp(double x) +{ + double res, resi, powx, e, factorial; + int i; + + if (isnan(x) != 0) { + return NAN; + } + + /* Values of x greater than 709.79 will cause overflow, returning INFINITY */ + if (x > 709.79) { + errno = ERANGE; + return HUGE_VAL; + } + + /* Get floor of exponent */ + x = modf(x, &e); + + /* Calculate most of the result */ + resi = quickPow(M_E, (int)e); + + /* Calculate rest of the result using Maclaurin series */ + factorial = 1.0; + powx = x; + res = 1.0; + + for (i = 2; i < 13; ++i) { + if (powx == 0.0) { + break; + } + res += powx / factorial; + factorial *= i; + powx *= x; + } + + return (res * resi); +} diff --git a/libm/phoenix/mathd/fabsd.c b/libm/phoenix/mathd/fabsd.c new file mode 100644 index 00000000..0c2d9bf8 --- /dev/null +++ b/libm/phoenix/mathd/fabsd.c @@ -0,0 +1,30 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include +#include "common.h" + + +double fabs(double x) +{ + if (isnan(x) != 0) { + return NAN; + } + + conv_t *conv = (conv_t *)&x; + conv->i.sign = 0; + + return x; +} diff --git a/libm/phoenix/mathd/floord.c b/libm/phoenix/mathd/floord.c new file mode 100644 index 00000000..77b9fed1 --- /dev/null +++ b/libm/phoenix/mathd/floord.c @@ -0,0 +1,34 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +double floor(double x) +{ + double ipart, fpart; + + if (isnan(x) != 0) { + return NAN; + } + + fpart = modf(x, &ipart); + + if ((x < 0.0) && ((fpart + x) != x)) { + ipart -= 1.0; + } + + return ipart; +} diff --git a/libm/phoenix/mathd/fmodd.c b/libm/phoenix/mathd/fmodd.c new file mode 100644 index 00000000..66059920 --- /dev/null +++ b/libm/phoenix/mathd/fmodd.c @@ -0,0 +1,42 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include +#include + + +double fmod(double number, double denom) +{ + double result, tquot; + + if (isnan(number) != 0 || isnan(denom) != 0) { + return NAN; + } + + if ((denom == 0.0) || (isinf(number) != 0)) { + errno = EDOM; + return NAN; + } + + if (((number == 0.0) && (denom != 0.0)) || + ((isinf(number) == 0) && (isinf(denom) != 0))) { + return number; + } + + modf(number / denom, &tquot); + result = tquot * denom; + + return number - result; +} diff --git a/libm/phoenix/mathd/frexpd.c b/libm/phoenix/mathd/frexpd.c new file mode 100644 index 00000000..88fb2d10 --- /dev/null +++ b/libm/phoenix/mathd/frexpd.c @@ -0,0 +1,45 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include +#include "common.h" + + +double frexp(double x, int *exp) +{ + if (isnan(x) != 0) { + return NAN; + } + + if (isinf(x) != 0) { + return x; + } + + if (x == 0.0) { + return x; + } + + conv_t *conv = (conv_t *)&x; + *exp = 0; + + if (conv->i.exponent == 0) { + normalizeSub(&x, exp); + } + + *exp += conv->i.exponent - 1022; + conv->i.exponent = 1022; + + return x; +} diff --git a/libm/phoenix/mathd/ldexpd.c b/libm/phoenix/mathd/ldexpd.c new file mode 100644 index 00000000..f9de327c --- /dev/null +++ b/libm/phoenix/mathd/ldexpd.c @@ -0,0 +1,55 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include +#include +#include "common.h" + + +double ldexp(double x, int exp) +{ + if (isnan(x) != 0) { + return NAN; + } + + if (x == 0.0) { + return x; + } + + conv_t *conv = (conv_t *)&x; + int exponent = 0; + + if (conv->i.exponent == 0) { + normalizeSub(&x, &exponent); + } + + exponent += conv->i.exponent + exp; + + if (exponent > 2046) { + errno = ERANGE; + return conv->i.sign ? -HUGE_VAL : HUGE_VAL; + } + + /* If result is subnormal */ + if (exponent < 0) { + createSub(&x, exponent); + conv->i.exponent = 0; + } + else { + conv->i.exponent = exponent; + } + + return x; +} diff --git a/libm/phoenix/mathd/log10d.c b/libm/phoenix/mathd/log10d.c new file mode 100644 index 00000000..b6119b6a --- /dev/null +++ b/libm/phoenix/mathd/log10d.c @@ -0,0 +1,23 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +/* Uses log10(x) = ln(x) / ln(10) identity */ +double log10(double x) +{ + return (log(x) / M_LN10); +} diff --git a/libm/phoenix/mathd/log2d.c b/libm/phoenix/mathd/log2d.c new file mode 100644 index 00000000..6905f55e --- /dev/null +++ b/libm/phoenix/mathd/log2d.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +double log2(double x) +{ + return (log(x) / M_LN2); +} diff --git a/libm/phoenix/mathd/logd.c b/libm/phoenix/mathd/logd.c new file mode 100644 index 00000000..1fc2d231 --- /dev/null +++ b/libm/phoenix/mathd/logd.c @@ -0,0 +1,63 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include +#include +#include "common.h" + + +double log(double x) +{ + double tmp, pow, res; + conv_t *conv = (conv_t *)&tmp; + int exp = 0, i; + + if (isnan(x) != 0) { + return NAN; + } + else if (x < 0.0) { + errno = EDOM; + return NAN; + } + else if (x == 0.0) { + errno = ERANGE; + return -HUGE_VAL; + } + else if (x == 1.0) { + return 0.0; + } + else if (isinf(x) != 0) { + return x; + } + + tmp = x; + + exp = conv->i.exponent - 1022; + + if (conv->i.exponent == 0) { + normalizeSub(&tmp, &exp); + } + + conv->i.exponent = 1022; + + tmp = (tmp - 1.0) / (tmp + 1.0); + + for (i = 1, res = 0.0, pow = tmp * tmp; i < 16; ++i) { + res += tmp / ((2 * i) - 1); + tmp *= pow; + } + + return ((2.0 * res) + (exp / M_LOG2E)); +} diff --git a/libm/phoenix/mathd/modfd.c b/libm/phoenix/mathd/modfd.c new file mode 100644 index 00000000..803b6082 --- /dev/null +++ b/libm/phoenix/mathd/modfd.c @@ -0,0 +1,57 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include +#include "common.h" + +double modf(double x, double *intpart) +{ + conv_t *conv = (conv_t *)&x; + double tmp = x; + int exp = conv->i.exponent - 1023; + uint64_t m, mask = 0xfffffffffffffLL; + + if (isnan(x) != 0) { + *intpart = NAN; + return NAN; + } + + if (exp > 52) { + *intpart = x; + return (conv->i.sign ? -0.0 : 0.0); + } + else if (exp < 0) { + *intpart = conv->i.sign ? -0.0 : 0.0; + return x; + } + + conv->i.mantisa = conv->i.mantisa & ~(mask >> exp); + *intpart = x; + x = tmp; + + m = conv->i.mantisa; + m &= mask >> exp; + + if (m == 0u) { + return 0.0; + } + + conv->i.mantisa = m & mask; + normalizeSub(&x, &exp); + + conv->i.exponent = exp + 1023; + + return x; +} diff --git a/math/power.c b/libm/phoenix/mathd/powd.c similarity index 56% rename from math/power.c rename to libm/phoenix/mathd/powd.c index 45837430..d9832e5e 100644 --- a/math/power.c +++ b/libm/phoenix/mathd/powd.c @@ -10,15 +10,13 @@ * * This file is part of Phoenix-RTOS. * - * %LICENSE% + * SPDX-License-Identifier: BSD-3-Clause */ -#include /* Needed for __ieee754_sqrt */ #include -#include -#include #include -#include +#include +#include #include "common.h" @@ -91,81 +89,3 @@ double pow(double base, double exponent) return res; } - - -double sqrt(double x) -{ - if (isnan(x) != 0) { - return NAN; - } - - if (x < 0.0) { - errno = EDOM; - return NAN; - } - - if ((x == 0.0) || (x == INFINITY)) { - return x; - } - -#ifdef __IEEE754_SQRT - return __ieee754_sqrt(x); -#else - /* Use reciprocal square root method: */ - double xn = 1.0 / x; - - /* IEEE-754 compliant: */ - conv_t *init = (conv_t *)&xn; - - /* +Infinity : */ - if ((init->i.mantisa == 0) && (init->i.exponent == 0x7FF)) { - return x; - } - - /* Subnormals: */ - if (init->i.exponent == 0) { - init->i.exponent = 0x1; - } - - /* Initial guess: */ - init->i.mantisa = (init->i.mantisa >> 1); - - if (init->i.exponent & 0x1) { - init->i.exponent = (init->i.exponent >> 1) + 0x200; - } - else { - init->i.exponent = (init->i.exponent >> 1) + 0x1FF; - init->i.mantisa |= (0x1ull << 51); - } - - /* Reciprocal sqare root iters (avoiding division): */ - for (int i = 0; i < 4; ++i) { - xn = xn * (1.5 - (0.5 * x * xn * xn)); - } - - return (xn * x); -#endif -} - - -float sqrtf(float x) -{ -#ifdef __IEEE754_SQRTF - if (isnan(x) != 0) { - return NAN; - } - - if (x < 0.0f) { - errno = EDOM; - return NAN; - } - - if ((x == 0.0f) || (x == INFINITY)) { - return x; - } - - return __ieee754_sqrtf(x); -#else - return (float)sqrt(x); -#endif -} diff --git a/libm/phoenix/mathd/roundd.c b/libm/phoenix/mathd/roundd.c new file mode 100644 index 00000000..7015edb3 --- /dev/null +++ b/libm/phoenix/mathd/roundd.c @@ -0,0 +1,37 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +double round(double x) +{ + double ret, frac; + + if (isnan(x) != 0) { + return NAN; + } + + frac = modf(x, &ret); + + if (frac >= 0.5) { + ret += 1.0; + } + else if (frac <= -0.5) { + ret -= 1.0; + } + + return ret; +} diff --git a/libm/phoenix/mathd/sind.c b/libm/phoenix/mathd/sind.c new file mode 100644 index 00000000..2f701521 --- /dev/null +++ b/libm/phoenix/mathd/sind.c @@ -0,0 +1,70 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * cos, sin, tan, acos, asin, atan + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Aleksander Kaminski, Jakub Smolaga + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include +#include + +/* Calculates value of sine using Maclaurin series via Horner's method. */ +double sin(double x) +{ + double res, xpow; + + if (isnan(x) != 0) { + return NAN; + } + + if (isinf(x) != 0) { + errno = EDOM; + return NAN; + } + + if (x == 0.0) { + return x; + } + + /* Normalize argument to -2*PI < x < 2*PI */ + x = fmod(x, 2.0 * M_PI); + + /* Normalize further to -PI < x < PI */ + if (x > M_PI) { + x -= 2 * M_PI; + } + else if (x < -M_PI) { + x += 2 * M_PI; + } + + if (x > M_PI_2) { + return cos(x - M_PI_2); + } + else if (x < -M_PI_2) { + return -cos(x + M_PI_2); + } + + xpow = x * x; + + res = 0x1.71b8ef6dcf572p-66; /* 1/21! */ + res = res * xpow - 0x1.2f49b46814157p-57; /* 1/19! */ + res = res * xpow + 0x1.952c77030ad4ap-49; /* 1/17! */ + res = res * xpow - 0x1.ae7f3e733b81fp-41; /* 1/15! */ + res = res * xpow + 0x1.6124613a86d09p-33; /* 1/13! */ + res = res * xpow - 0x1.ae64567f544e4p-26; /* 1/11! */ + res = res * xpow + 0x1.71de3a556c734p-19; /* 1/9! */ + res = res * xpow - 0x1.a01a01a01a01ap-13; /* 1/7! */ + res = res * xpow + 0x1.1111111111111p-7; /* 1/5! */ + res = res * xpow - 0x1.5555555555555p-3; /* 1/3! */ + res = res * xpow + 0x1.0p0; /* 1/1! */ + + return res * x; +} diff --git a/math/hyper.c b/libm/phoenix/mathd/sinhd.c similarity index 60% rename from math/hyper.c rename to libm/phoenix/mathd/sinhd.c index c3afb86f..c17e310b 100644 --- a/math/hyper.c +++ b/libm/phoenix/mathd/sinhd.c @@ -10,35 +10,12 @@ * * This file is part of Phoenix-RTOS. * - * %LICENSE% + * SPDX-License-Identifier: BSD-3-Clause */ - + #include -#include #include - - -double cosh(double x) -{ - double y; - - if (isnan(x) != 0) { - return NAN; - } - - if (isinf(x) != 0) { - return INFINITY; - } - - /* Make sure exp(x) is not infinity */ - if (x < 709.78) { - return ((exp(x) + exp(-x)) / 2.0); - } - else { - y = cosh(x / 2.0); - return ((2.0 * y * y) - 1.0); - } -} +#include double sinh(double x) @@ -70,21 +47,3 @@ double sinh(double x) return f; } } - - -double tanh(double x) -{ - if (isnan(x) != 0) { - return NAN; - } - - if (x == INFINITY) { - return 1.0; - } - else if (x == -INFINITY) { - return -1.0; - } - - /* cosh is never equal to zero */ - return (sinh(x) / cosh(x)); -} diff --git a/libm/phoenix/mathd/sqrtd.c b/libm/phoenix/mathd/sqrtd.c new file mode 100644 index 00000000..50425050 --- /dev/null +++ b/libm/phoenix/mathd/sqrtd.c @@ -0,0 +1,69 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * pow, sqrt + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include +#include +#include "common.h" + + +double sqrt(double x) +{ + if (isnan(x) != 0) { + return NAN; + } + + if (x < 0.0) { + errno = EDOM; + return NAN; + } + + if ((x == 0.0) || (x == INFINITY)) { + return x; + } + + /* Use reciprocal square root method: */ + double xn = 1.0 / x; + + /* IEEE-754 compliant: */ + conv_t *init = (conv_t *)&xn; + + /* +Infinity : */ + if ((init->i.mantisa == 0) && (init->i.exponent == 0x7FF)) { + return x; + } + + /* Subnormals: */ + if (init->i.exponent == 0) { + init->i.exponent = 0x1; + } + + /* Initial guess: */ + init->i.mantisa = (init->i.mantisa >> 1); + + if (init->i.exponent & 0x1) { + init->i.exponent = (init->i.exponent >> 1) + 0x200; + } + else { + init->i.exponent = (init->i.exponent >> 1) + 0x1FF; + init->i.mantisa |= (0x1ull << 51); + } + + /* Reciprocal sqare root iters (avoiding division): */ + for (int i = 0; i < 4; ++i) { + xn = xn * (1.5 - (0.5 * x * xn * xn)); + } + + return (xn * x); +} diff --git a/libm/phoenix/mathd/tand.c b/libm/phoenix/mathd/tand.c new file mode 100644 index 00000000..646297a6 --- /dev/null +++ b/libm/phoenix/mathd/tand.c @@ -0,0 +1,41 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * cos, sin, tan, acos, asin, atan + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Aleksander Kaminski, Jakub Smolaga + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include +#include + + +double tan(double x) +{ + double c; + + if (isnan(x) != 0) { + return NAN; + } + + if (isinf(x) != 0) { + errno = EDOM; + return NAN; + } + + c = cos(x); + + if (c != 0.0) { + return (sin(x) / c); + } + else { + return 0.0; + } +} diff --git a/libm/phoenix/mathd/tanhd.c b/libm/phoenix/mathd/tanhd.c new file mode 100644 index 00000000..dbb840c6 --- /dev/null +++ b/libm/phoenix/mathd/tanhd.c @@ -0,0 +1,34 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * hyberbolic functions + * + * Copyright 2018 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +double tanh(double x) +{ + if (isnan(x) != 0) { + return NAN; + } + + if (x == INFINITY) { + return 1.0; + } + else if (x == -INFINITY) { + return -1.0; + } + + /* cosh is never equal to zero */ + return (sinh(x) / cosh(x)); +} diff --git a/libm/phoenix/mathd/truncd.c b/libm/phoenix/mathd/truncd.c new file mode 100644 index 00000000..b5678d3b --- /dev/null +++ b/libm/phoenix/mathd/truncd.c @@ -0,0 +1,30 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +double trunc(double x) +{ + double ret; + + if (isnan(x) != 0) { + return NAN; + } + + modf(x, &ret); + + return ret; +} diff --git a/libm/phoenix/mathf/acosf.c b/libm/phoenix/mathf/acosf.c new file mode 100644 index 00000000..d90567be --- /dev/null +++ b/libm/phoenix/mathf/acosf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * cos, sin, tan, acos, asin, atan + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Aleksander Kaminski, Jakub Smolaga + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float acosf(float x) +{ + return (float)acos((double)x); +} diff --git a/libm/phoenix/mathf/asinf.c b/libm/phoenix/mathf/asinf.c new file mode 100644 index 00000000..78c140bc --- /dev/null +++ b/libm/phoenix/mathf/asinf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * cos, sin, tan, acos, asin, atan + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Aleksander Kaminski, Jakub Smolaga + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float asinf(float x) +{ + return (float)asin((double)x); +} diff --git a/libm/phoenix/mathf/atan2f.c b/libm/phoenix/mathf/atan2f.c new file mode 100644 index 00000000..d19ae29a --- /dev/null +++ b/libm/phoenix/mathf/atan2f.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * cos, sin, tan, acos, asin, atan + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Aleksander Kaminski, Jakub Smolaga + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float atan2f(float y, float x) +{ + return (float)atan2((double)y, (double)x); +} diff --git a/libm/phoenix/mathf/atanf.c b/libm/phoenix/mathf/atanf.c new file mode 100644 index 00000000..852b75c0 --- /dev/null +++ b/libm/phoenix/mathf/atanf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * cos, sin, tan, acos, asin, atan + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Aleksander Kaminski, Jakub Smolaga + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float atanf(float x) +{ + return (float)atan((double)x); +} diff --git a/libm/phoenix/mathf/ceilf.c b/libm/phoenix/mathf/ceilf.c new file mode 100644 index 00000000..c600a9f7 --- /dev/null +++ b/libm/phoenix/mathf/ceilf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float ceilf(float x) +{ + return (float)ceil(x); +} diff --git a/libm/phoenix/mathf/cosf.c b/libm/phoenix/mathf/cosf.c new file mode 100644 index 00000000..e3ae08a7 --- /dev/null +++ b/libm/phoenix/mathf/cosf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * cos, sin, tan, acos, asin, atan + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Aleksander Kaminski, Jakub Smolaga + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float cosf(float x) +{ + return (float)cos((double)x); +} diff --git a/libm/phoenix/mathf/coshf.c b/libm/phoenix/mathf/coshf.c new file mode 100644 index 00000000..35097084 --- /dev/null +++ b/libm/phoenix/mathf/coshf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * hyberbolic functions + * + * Copyright 2018 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float coshf(float x) +{ + return (float)cosh((double)x); +} diff --git a/libm/phoenix/mathf/expf.c b/libm/phoenix/mathf/expf.c new file mode 100644 index 00000000..51925724 --- /dev/null +++ b/libm/phoenix/mathf/expf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float expf(float x) +{ + return (float)exp((double)x); +} diff --git a/libm/phoenix/mathf/fabsf.c b/libm/phoenix/mathf/fabsf.c new file mode 100644 index 00000000..290fc2ed --- /dev/null +++ b/libm/phoenix/mathf/fabsf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float fabsf(float x) +{ + return (float)fabs(x); +} diff --git a/libm/phoenix/mathf/floorf.c b/libm/phoenix/mathf/floorf.c new file mode 100644 index 00000000..948bd041 --- /dev/null +++ b/libm/phoenix/mathf/floorf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float floorf(float x) +{ + return (float)floor(x); +} diff --git a/libm/phoenix/mathf/fmodf.c b/libm/phoenix/mathf/fmodf.c new file mode 100644 index 00000000..1a750ec1 --- /dev/null +++ b/libm/phoenix/mathf/fmodf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float fmodf(float x, float y) +{ + return (float)fmod((double)x, (double)y); +} diff --git a/libm/phoenix/mathf/frexpf.c b/libm/phoenix/mathf/frexpf.c new file mode 100644 index 00000000..f6d45403 --- /dev/null +++ b/libm/phoenix/mathf/frexpf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float frexpf(float x, int *exp) +{ + return (float)frexp((double)x, exp); +} diff --git a/libm/phoenix/mathf/ldexpf.c b/libm/phoenix/mathf/ldexpf.c new file mode 100644 index 00000000..246ce5b6 --- /dev/null +++ b/libm/phoenix/mathf/ldexpf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float ldexpf(float x, int exp) +{ + return (float)ldexp((double)x, exp); +} diff --git a/libm/phoenix/mathf/log10f.c b/libm/phoenix/mathf/log10f.c new file mode 100644 index 00000000..0e5de3b1 --- /dev/null +++ b/libm/phoenix/mathf/log10f.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float log10f(float x) +{ + return (float)log10((double)x); +} diff --git a/libm/phoenix/mathf/logf.c b/libm/phoenix/mathf/logf.c new file mode 100644 index 00000000..cc8f282b --- /dev/null +++ b/libm/phoenix/mathf/logf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float logf(float x) +{ + return (float)log((double)x); +} diff --git a/libm/phoenix/mathf/modff.c b/libm/phoenix/mathf/modff.c new file mode 100644 index 00000000..162a4e3a --- /dev/null +++ b/libm/phoenix/mathf/modff.c @@ -0,0 +1,27 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float modff(float x, float *intpart) +{ + double ret, tmp; + + ret = modf(x, &tmp); + *intpart = tmp; + + return ret; +} diff --git a/libm/phoenix/mathf/powf.c b/libm/phoenix/mathf/powf.c new file mode 100644 index 00000000..985ad80f --- /dev/null +++ b/libm/phoenix/mathf/powf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * pow, sqrt + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float powf(float x, float exp) +{ + return (float)pow((double)x, (double)exp); +} diff --git a/libm/phoenix/mathf/roundf.c b/libm/phoenix/mathf/roundf.c new file mode 100644 index 00000000..b91a3d29 --- /dev/null +++ b/libm/phoenix/mathf/roundf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float roundf(float x) +{ + return (float)round(x); +} diff --git a/libm/phoenix/mathf/sinf.c b/libm/phoenix/mathf/sinf.c new file mode 100644 index 00000000..864fa2ed --- /dev/null +++ b/libm/phoenix/mathf/sinf.c @@ -0,0 +1,21 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * cos, sin, tan, acos, asin, atan + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Aleksander Kaminski, Jakub Smolaga + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + +float sinf(float x) +{ + return (float)sin((double)x); +} diff --git a/libm/phoenix/mathf/sinhf.c b/libm/phoenix/mathf/sinhf.c new file mode 100644 index 00000000..739312b2 --- /dev/null +++ b/libm/phoenix/mathf/sinhf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * hyberbolic functions + * + * Copyright 2018 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float sinhf(float x) +{ + return (double)sinh((double)x); +} diff --git a/libm/phoenix/mathf/sqrtf.c b/libm/phoenix/mathf/sqrtf.c new file mode 100644 index 00000000..660832eb --- /dev/null +++ b/libm/phoenix/mathf/sqrtf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * pow, sqrt + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float sqrtf(float x) +{ + return (float)sqrt(x); +} diff --git a/libm/phoenix/mathf/tanf.c b/libm/phoenix/mathf/tanf.c new file mode 100644 index 00000000..bfe46bdc --- /dev/null +++ b/libm/phoenix/mathf/tanf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * cos, sin, tan, acos, asin, atan + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Aleksander Kaminski, Jakub Smolaga + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float tanf(float x) +{ + return (float)tan((double)x); +} diff --git a/libm/phoenix/mathf/tanhf.c b/libm/phoenix/mathf/tanhf.c new file mode 100644 index 00000000..b9a49f35 --- /dev/null +++ b/libm/phoenix/mathf/tanhf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * hyberbolic functions + * + * Copyright 2018 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float tanhf(float x) +{ + return (float)tanh((double)x); +} diff --git a/libm/phoenix/mathf/truncf.c b/libm/phoenix/mathf/truncf.c new file mode 100644 index 00000000..d1200e3d --- /dev/null +++ b/libm/phoenix/mathf/truncf.c @@ -0,0 +1,22 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs + * + * Copyright 2017 Phoenix Systems + * Author: Aleksander Kaminski + * + * This file is part of Phoenix-RTOS. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include + + +float truncf(float x) +{ + return (float)trunc(x); +} diff --git a/math/Makefile b/math/Makefile deleted file mode 100644 index f8f1a967..00000000 --- a/math/Makefile +++ /dev/null @@ -1,7 +0,0 @@ -# -# Makefile for libphoenix/math -# -# Copyright 2012-2015, 2020 Phoenix Systems -# - -OBJS += $(addprefix $(PREFIX_O)math/, trig.o power.o exp.o common.o complex.o hyper.o) diff --git a/math/exp.c b/math/exp.c deleted file mode 100644 index 77c9c76a..00000000 --- a/math/exp.c +++ /dev/null @@ -1,411 +0,0 @@ -/* - * Phoenix-RTOS - * - * libphoenix - * - * exp, frexp, ldexp, log, log10, modf, ceil, floor, fmod, fabs - * - * Copyright 2017 Phoenix Systems - * Author: Aleksander Kaminski - * - * This file is part of Phoenix-RTOS. - * - * %LICENSE% - */ - -#include -#include -#include -#include "common.h" - - -double frexp(double x, int *exp) -{ - if (isnan(x) != 0) { - return NAN; - } - - if (isinf(x) != 0) { - return x; - } - - if (x == 0.0) { - return x; - } - - conv_t *conv = (conv_t *)&x; - *exp = 0; - - if (conv->i.exponent == 0) { - normalizeSub(&x, exp); - } - - *exp += conv->i.exponent - 1022; - conv->i.exponent = 1022; - - return x; -} - - -double ldexp(double x, int exp) -{ - if (isnan(x) != 0) { - return NAN; - } - - if (x == 0.0) { - return x; - } - - conv_t *conv = (conv_t *)&x; - int exponent = 0; - - if (conv->i.exponent == 0) { - normalizeSub(&x, &exponent); - } - - exponent += conv->i.exponent + exp; - - if (exponent > 2046) { - errno = ERANGE; - return conv->i.sign ? -HUGE_VAL : HUGE_VAL; - } - - /* If result is subnormal */ - if (exponent < 0) { - createSub(&x, exponent); - conv->i.exponent = 0; - } - else { - conv->i.exponent = exponent; - } - - return x; -} - - -double log(double x) -{ - double tmp, pow, res; - conv_t *conv = (conv_t *)&tmp; - int exp = 0, i; - - if (isnan(x) != 0) { - return NAN; - } - else if (x < 0.0) { - errno = EDOM; - return NAN; - } - else if (x == 0.0) { - errno = ERANGE; - return -HUGE_VAL; - } - else if (x == 1.0) { - return 0.0; - } - else if (isinf(x) != 0) { - return x; - } - - tmp = x; - - exp = conv->i.exponent - 1022; - - if (conv->i.exponent == 0) { - normalizeSub(&tmp, &exp); - } - - conv->i.exponent = 1022; - - tmp = (tmp - 1.0) / (tmp + 1.0); - - for (i = 1, res = 0.0, pow = tmp * tmp; i < 16; ++i) { - res += tmp / ((2 * i) - 1); - tmp *= pow; - } - - return ((2.0 * res) + (exp / M_LOG2E)); -} - - -double log2(double x) -{ - return (log(x) / M_LN2); -} - - -/* Uses log10(x) = ln(x) / ln(10) identity */ -double log10(double x) -{ - return (log(x) / M_LN10); -} - - -double modf(double x, double *intpart) -{ - conv_t *conv = (conv_t *)&x; - double tmp = x; - int exp = conv->i.exponent - 1023; - uint64_t m, mask = 0xfffffffffffffLL; - - if (isnan(x) != 0) { - *intpart = NAN; - return NAN; - } - - if (exp > 52) { - *intpart = x; - return (conv->i.sign ? -0.0 : 0.0); - } - else if (exp < 0) { - *intpart = conv->i.sign ? -0.0 : 0.0; - return x; - } - - conv->i.mantisa = conv->i.mantisa & ~(mask >> exp); - *intpart = x; - x = tmp; - - m = conv->i.mantisa; - m &= mask >> exp; - - if (m == 0u) { - return 0.0; - } - - conv->i.mantisa = m & mask; - normalizeSub(&x, &exp); - - conv->i.exponent = exp + 1023; - - return x; -} - -float modff(float x, float *intpart) -{ - double ret, tmp; - - ret = modf(x, &tmp); - *intpart = tmp; - - return ret; -} - -/* Uses quick powering and Maclaurin series to calculate value of e^x */ -double exp(double x) -{ - double res, resi, powx, e, factorial; - int i; - - if (isnan(x) != 0) { - return NAN; - } - - /* Values of x greater than 709.79 will cause overflow, returning INFINITY */ - if (x > 709.79) { - errno = ERANGE; - return HUGE_VAL; - } - - /* Get floor of exponent */ - x = modf(x, &e); - - /* Calculate most of the result */ - resi = quickPow(M_E, (int)e); - - /* Calculate rest of the result using Maclaurin series */ - factorial = 1.0; - powx = x; - res = 1.0; - - for (i = 2; i < 13; ++i) { - if (powx == 0.0) { - break; - } - res += powx / factorial; - factorial *= i; - powx *= x; - } - - return (res * resi); -} - - -double ceil(double x) -{ -#ifdef __IEEE754_CEIL - return __ieee754_ceil(x); -#else - double ipart, fpart; - - if (isnan(x) != 0) { - return NAN; - } - - fpart = modf(x, &ipart); - - if ((x > 0.0) && ((fpart + x) != x)) { - ipart += 1.0; - } - - return ipart; -#endif -} - - -float ceilf(float x) -{ -#ifdef __IEEE754_CEILF - return __ieee754_ceilf(x); -#else - return (float)ceil(x); -#endif -} - - -double floor(double x) -{ -#ifdef __IEEE754_FLOOR - return __ieee754_floor(x); -#else - double ipart, fpart; - - if (isnan(x) != 0) { - return NAN; - } - - fpart = modf(x, &ipart); - - if ((x < 0.0) && ((fpart + x) != x)) { - ipart -= 1.0; - } - - return ipart; -#endif -} - - -float floorf(float x) -{ -#ifdef __IEEE754_FLOORF - return __ieee754_floorf(x); -#else - return (float)floor(x); -#endif -} - - -double fmod(double number, double denom) -{ - double result, tquot; - - if (isnan(number) != 0 || isnan(denom) != 0) { - return NAN; - } - - if ((denom == 0.0) || (isinf(number) != 0)) { - errno = EDOM; - return NAN; - } - - if (((number == 0.0) && (denom != 0.0)) || - ((isinf(number) == 0) && (isinf(denom) != 0))) { - return number; - } - - modf(number / denom, &tquot); - result = tquot * denom; - - return number - result; -} - - -double round(double x) -{ -#ifdef __IEEE754_ROUND - return __ieee754_round(x); -#else - double ret, frac; - - if (isnan(x) != 0) { - return NAN; - } - - frac = modf(x, &ret); - - if (frac >= 0.5) { - ret += 1.0; - } - else if (frac <= -0.5) { - ret -= 1.0; - } - - return ret; -#endif -} - - -float roundf(float x) -{ -#ifdef __IEEE754_ROUNDF - return __ieee754_roundf(x); -#else - return (float)round(x); -#endif -} - - -double trunc(double x) -{ -#ifdef __IEEE754_TRUNC - return __ieee754_trunc(x); -#else - double ret; - - if (isnan(x) != 0) { - return NAN; - } - - modf(x, &ret); - - return ret; -#endif -} - - -float truncf(float x) -{ -#ifdef __IEEE754_TRUNCF - return __ieee754_truncf(x); -#else - return (float)trunc(x); -#endif -} - - -double fabs(double x) -{ -#ifdef __IEEE754_FABS - return __ieee754_fabs(x); -#else - if (isnan(x) != 0) { - return NAN; - } - - conv_t *conv = (conv_t *)&x; - conv->i.sign = 0; - - return x; -#endif -} - - -float fabsf(float x) -{ -#ifdef __IEEE754_FABSF - return __ieee754_fabsf(x); -#else - return (float)fabs(x); -#endif -} diff --git a/math/trig.c b/math/trig.c deleted file mode 100644 index ac1008e9..00000000 --- a/math/trig.c +++ /dev/null @@ -1,293 +0,0 @@ -/* - * Phoenix-RTOS - * - * libphoenix - * - * cos, sin, tan, acos, asin, atan - * - * Copyright 2017, 2018, 2026 Phoenix Systems - * Author: Aleksander Kaminski, Jakub Smolaga - * - * This file is part of Phoenix-RTOS. - * - * %LICENSE% - */ - -#include -#include - - -/* Calculates value of cosine using Maclaurin series via Horner's method. */ -double cos(double x) -{ - double res, xpow; - - if (isnan(x) != 0) { - return NAN; - } - - if (isinf(x) != 0) { - errno = EDOM; - return NAN; - } - - /* Normalize argument to -2*PI < x < 2*PI */ - x = fmod(x, 2.0 * M_PI); - - /* Normalize further to -PI < x < PI */ - if (x > M_PI) { - x -= 2.0 * M_PI; - } - else if (x < -M_PI) { - x += 2.0 * M_PI; - } - - if (x > M_PI_2) { - return -sin(x - M_PI_2); - } - else if (x < -M_PI_2) { - return sin(x + M_PI_2); - } - - xpow = x * x; - - res = 0x1.e542ba4020225p-62; /* 1/20! */ - res = res * xpow - 0x1.6827863b97d97p-53; /* 1/18! */ - res = res * xpow + 0x1.ae7f3e733b81fp-45; /* 1/16! */ - res = res * xpow - 0x1.93974a8c07c9dp-37; /* 1/14! */ - res = res * xpow + 0x1.1eed8eff8d898p-29; /* 1/12! */ - res = res * xpow - 0x1.27e4fb7789f5cp-22; /* 1/10! */ - res = res * xpow + 0x1.a01a01a01a01ap-16; /* 1/8! */ - res = res * xpow - 0x1.6c16c16c16c17p-10; /* 1/6! */ - res = res * xpow + 0x1.5555555555555p-5; /* 1/4! */ - res = res * xpow - 0x1p-1; /* 1/2! */ - res = res * xpow + 0x1p0; /* 1/0! */ - - return res; -} - - -/* Calculates value of sine using Maclaurin series via Horner's method. */ -double sin(double x) -{ - double res, xpow; - - if (isnan(x) != 0) { - return NAN; - } - - if (isinf(x) != 0) { - errno = EDOM; - return NAN; - } - - if (x == 0.0) { - return x; - } - - /* Normalize argument to -2*PI < x < 2*PI */ - x = fmod(x, 2.0 * M_PI); - - /* Normalize further to -PI < x < PI */ - if (x > M_PI) { - x -= 2 * M_PI; - } - else if (x < -M_PI) { - x += 2 * M_PI; - } - - if (x > M_PI_2) { - return cos(x - M_PI_2); - } - else if (x < -M_PI_2) { - return -cos(x + M_PI_2); - } - - xpow = x * x; - - res = 0x1.71b8ef6dcf572p-66; /* 1/21! */ - res = res * xpow - 0x1.2f49b46814157p-57; /* 1/19! */ - res = res * xpow + 0x1.952c77030ad4ap-49; /* 1/17! */ - res = res * xpow - 0x1.ae7f3e733b81fp-41; /* 1/15! */ - res = res * xpow + 0x1.6124613a86d09p-33; /* 1/13! */ - res = res * xpow - 0x1.ae64567f544e4p-26; /* 1/11! */ - res = res * xpow + 0x1.71de3a556c734p-19; /* 1/9! */ - res = res * xpow - 0x1.a01a01a01a01ap-13; /* 1/7! */ - res = res * xpow + 0x1.1111111111111p-7; /* 1/5! */ - res = res * xpow - 0x1.5555555555555p-3; /* 1/3! */ - res = res * xpow + 0x1.0p0; /* 1/1! */ - - return res * x; -} - - -double tan(double x) -{ - double c; - - if (isnan(x) != 0) { - return NAN; - } - - if (isinf(x) != 0) { - errno = EDOM; - return NAN; - } - - c = cos(x); - - if (c != 0.0) { - return (sin(x) / c); - } - else { - return 0.0; - } -} - - -/* Calculates value of arc cosine using secant method */ -double acos(double x) -{ - double xa = 0, xb = M_PI, ya, yb, t; - int i; - - if (isnan(x) != 0) { - return NAN; - } - - if ((x > 1.0) || (x < -1.0)) { - errno = EDOM; - return NAN; - } - - for (i = 0; i < 16; ++i) { - ya = cos(xa) - x; - yb = cos(xb) - x; - - t = ya - yb; - - if (t == 0.0) { - break; - } - - t = ((ya * xb) - (yb * xa)) / t; - xa = xb; - xb = t; - } - - return xb; -} - - -/* Calculates value of arc sine using asin(x) = pi/2 - acos(x) relationship. */ -double asin(double x) -{ - if (isnan(x) != 0) { - return NAN; - } - - if (x == 0.0) { - return x; - } - - return M_PI_2 - acos(x); -} - - -static const double atan_wi[] = { 0.2152638534631578, 0.2152638534631578, 0.2051984637212956, - 0.2051984637212956, 0.1855383974779378, 0.1855383974779378, 0.1572031671581935, - 0.1572031671581935, 0.1215185706879032, 0.1215185706879032, 0.0801580871597602, - 0.0801580871597602, 0.0351194603317519, 0.0351194603317519 }; - - -static const double atan_xi[] = { -0.1080549487073437, 0.1080549487073437, -0.3191123689278897, - 0.3191123689278897, -0.5152486363581541, 0.5152486363581541, -0.6872929048116855, - 0.6872929048116855, -0.8272013150697650, 0.8272013150697650, -0.9284348836635735, - 0.9284348836635735, -0.9862838086968123, 0.9862838086968123 }; - - -double atan(double x) -{ - double res = 1.0, h, a; - int i, s; - - if (isnan(x) != 0) { - return NAN; - } - - s = (x < 0.0) ? -1 : 1; - x *= (double)s; - h = x / 2; - - if (x == 0.0) { - return x; - } - - if (x > 1.0) { - return ((M_PI_2 - atan(1.0 / x)) * (double)s); - } - - for (i = 0, res = 0.0; i < sizeof(atan_wi) / sizeof(atan_wi[0]); ++i) { - a = (h * atan_xi[i]) + h; - res += atan_wi[i] / ((a * a) + 1); - } - - return (res * h * (double)s); -} - - -double atan2(double y, double x) -{ - if ((isnan(y) != 0) || (isnan(x) != 0)) { - return NAN; - } - - if (y == 0.0) { - if (x < 0.0 || signbit(x) != 0) { - // x is negative or -0.0 - return (signbit(y) != 0) ? -M_PI : M_PI; - } - else { - // x is positive or +0.0 - return y; - } - } - - if ((isinf(y) != 0) && (isinf(x) != 0)) { - int sy = signbit(y); - int sx = signbit(x); - - if ((sy == 0) && (sx != 0)) { - return M_PI_4 * 3; - } - else if ((sy != 0) && (sx != 0)) { - return -M_PI_4 * 3; - } - else if ((sy == 0) && (sx == 0)) { - return M_PI_4; - } - else if ((sy != 0) && (sx == 0)) { - return -M_PI_4; - } - } - - if (x > 0.0) { - return atan(y / x); - } - else if (x < 0.0) { - if (y >= 0.0) { - return atan(y / x) + M_PI; - } - - return atan(y / x) - M_PI; - } - - if (y > 0.0) { - return M_PI_2; - } - else if (y < 0.0) { - return -M_PI_2; - } - - return 0.0; -} From 0473f494c71d567219338e706711ec507dcf2a96 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miko=C5=82aj=20Mata=C5=82owski?= Date: Fri, 17 Jul 2026 17:28:20 +0200 Subject: [PATCH 2/3] libm: Add libmcs implementation YT: RTOS-1132 --- .gitmodules | 4 ++++ libm/README.md | 38 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 42 insertions(+) create mode 100644 .gitmodules create mode 100644 libm/README.md diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 00000000..cf947174 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,4 @@ +[submodule "libm/libmcs"] + path = libm/libmcs + url = https://github.com/phoenix-rtos/libmcs.git + branch = phoenix diff --git a/libm/README.md b/libm/README.md new file mode 100644 index 00000000..82f41611 --- /dev/null +++ b/libm/README.md @@ -0,0 +1,38 @@ +# Libm +This directory contains implementation of math library. It currently links to libphoenix with intention to be compiled to separate libm.a library. There are two implementations for the following cases: +- libphoenix implementation is not complete. +- version for critical systems is nice to have :) + +Usage: For each project it is suggested to set following +- ```LIBM_USE_LIBMCS ?= y/n``` - select implementation of math library for critical systems (if not set will use phoenix implementation) +- ```LIBM_WANT_COMPLEX ?= y/n``` - select if project wants implementation of complex routines +- ```LIBM_LIBMCS_DAZ ?= y/n``` - (applicable only for libmcs) select whether denormals are zero. This is especially useful if underlying FPU has limited functionalities implemented. +- ```LIBM_USE_HW ?= y``` - enable hardware intrinsics for math functions + +# Architecture specific code +Some targets might support instructions that implement parts of `math.h` required functions. Such implementation may be placed `arch/`, and they will take precedence before software implementation, unless stated otherwise via appropriate flags. Please refer to already existing code if new architecture specific implementation is to be added. + +# LIBMCS +This is a copy of libmcs math library implementation from https://gitlab.com/gtd-gmbh/libmcs (v1.3.1) + +For documentation please refer to /docs directory present on the main repo linked above. + +# Phoenix +This is implementation of math library developed alongside Pheonix RTOS operating system. Keep in mind that this implementation is not complete, not all functions in C99 standard are implemented. + +# Notes +1. This is a temporary state, moving math library to be separate to libphoenix +2. Some hardware instrinsics for math functions are still not implemented + + + + + + + + + + + + + From c9483367f0e305435e7064049d9c4b04c306aee4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miko=C5=82aj=20Mata=C5=82owski?= Date: Fri, 17 Jul 2026 17:28:48 +0200 Subject: [PATCH 3/3] libm: Add architecture specific math implementation YT: RTOS-1132 --- libm/arch/Makefile | 12 ++++++++ libm/arch/aarch64/ceild.c | 27 +++++++++++++++++ libm/arch/aarch64/ceilf.c | 27 +++++++++++++++++ libm/arch/aarch64/fabsd.c | 30 +++++++++++++++++++ libm/arch/aarch64/fabsf.c | 30 +++++++++++++++++++ libm/arch/aarch64/floord.c | 27 +++++++++++++++++ libm/arch/aarch64/floorf.c | 27 +++++++++++++++++ libm/arch/aarch64/roundd.c | 27 +++++++++++++++++ libm/arch/aarch64/roundf.c | 27 +++++++++++++++++ libm/arch/aarch64/sqrtd.c | 37 +++++++++++++++++++++++ libm/arch/aarch64/sqrtf.c | 37 +++++++++++++++++++++++ libm/arch/aarch64/truncd.c | 27 +++++++++++++++++ libm/arch/aarch64/truncf.c | 27 +++++++++++++++++ libm/arch/armv7a/sqrtd.c | 36 +++++++++++++++++++++++ libm/arch/armv7a/sqrtf.c | 36 +++++++++++++++++++++++ libm/arch/armv7m/.gitkeep | 0 libm/arch/armv7r/sqrtd.c | 36 +++++++++++++++++++++++ libm/arch/armv7r/sqrtf.c | 36 +++++++++++++++++++++++ libm/arch/armv8m/.gitkeep | 0 libm/arch/armv8r/sqrtd.c | 37 +++++++++++++++++++++++ libm/arch/armv8r/sqrtf.c | 37 +++++++++++++++++++++++ libm/arch/ia32/sqrtd.c | 54 ++++++++++++++++++++++++++++++++++ libm/arch/ia32/sqrtf.c | 36 +++++++++++++++++++++++ libm/arch/riscv64/sqrtd.c | 39 ++++++++++++++++++++++++ libm/arch/riscv64/sqrtf.c | 40 +++++++++++++++++++++++++ libm/arch/sparcv8leon/.gitkeep | 0 26 files changed, 749 insertions(+) create mode 100644 libm/arch/Makefile create mode 100644 libm/arch/aarch64/ceild.c create mode 100644 libm/arch/aarch64/ceilf.c create mode 100644 libm/arch/aarch64/fabsd.c create mode 100644 libm/arch/aarch64/fabsf.c create mode 100644 libm/arch/aarch64/floord.c create mode 100644 libm/arch/aarch64/floorf.c create mode 100644 libm/arch/aarch64/roundd.c create mode 100644 libm/arch/aarch64/roundf.c create mode 100644 libm/arch/aarch64/sqrtd.c create mode 100644 libm/arch/aarch64/sqrtf.c create mode 100644 libm/arch/aarch64/truncd.c create mode 100644 libm/arch/aarch64/truncf.c create mode 100644 libm/arch/armv7a/sqrtd.c create mode 100644 libm/arch/armv7a/sqrtf.c create mode 100644 libm/arch/armv7m/.gitkeep create mode 100644 libm/arch/armv7r/sqrtd.c create mode 100644 libm/arch/armv7r/sqrtf.c create mode 100644 libm/arch/armv8m/.gitkeep create mode 100644 libm/arch/armv8r/sqrtd.c create mode 100644 libm/arch/armv8r/sqrtf.c create mode 100644 libm/arch/ia32/sqrtd.c create mode 100644 libm/arch/ia32/sqrtf.c create mode 100644 libm/arch/riscv64/sqrtd.c create mode 100644 libm/arch/riscv64/sqrtf.c create mode 100644 libm/arch/sparcv8leon/.gitkeep diff --git a/libm/arch/Makefile b/libm/arch/Makefile new file mode 100644 index 00000000..dcbfeb4e --- /dev/null +++ b/libm/arch/Makefile @@ -0,0 +1,12 @@ +# +# Makefile for target specific math +# +# Copyright 2026 Phoenix Systems +# Author: Mikolaj Matalowski +# +# This file is part of Phoenix-RTOS. +# +# %LICENSE% +# + +HW_OBJS := $(addprefix $(PREFIX_O), $(patsubst %.c,%.o,$(wildcard libm/arch/$(TARGET_SUFF)/*.c))) \ No newline at end of file diff --git a/libm/arch/aarch64/ceild.c b/libm/arch/aarch64/ceild.c new file mode 100644 index 00000000..710bf1e8 --- /dev/null +++ b/libm/arch/aarch64/ceild.c @@ -0,0 +1,27 @@ +/* + * Phoenix-RTOS + * + * libphoenix/libm + * + * Architecture dependent part (arch/aarch64) + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Aleksander Kaminski, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ + +#include + + +double ceil(double x) +{ + if (isnan(x) != 0) { + return NAN; + } + + __asm__ volatile("frintp %d0, %d1" : "=w"(x) : "w"(x)); + return x; +} diff --git a/libm/arch/aarch64/ceilf.c b/libm/arch/aarch64/ceilf.c new file mode 100644 index 00000000..e17e4f77 --- /dev/null +++ b/libm/arch/aarch64/ceilf.c @@ -0,0 +1,27 @@ +/* + * Phoenix-RTOS + * + * libphoenix/libm + * + * Architecture dependent part (arch/aarch64) + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Aleksander Kaminski, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ + +#include + + +float ceilf(float x) +{ + if (isnan(x) != 0) { + return NAN; + } + + __asm__ volatile("frintp %s0, %s1" : "=w"(x) : "w"(x)); + return x; +} diff --git a/libm/arch/aarch64/fabsd.c b/libm/arch/aarch64/fabsd.c new file mode 100644 index 00000000..91bbb44f --- /dev/null +++ b/libm/arch/aarch64/fabsd.c @@ -0,0 +1,30 @@ +/* + * Phoenix-RTOS + * + * libphoenix/libm + * + * Architecture dependent part (arch/aarch64) + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Aleksander Kaminski, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ + +#include + + +double fabs(double x) +{ + if (isnan(x)) { + return x; + } + if (0 < x) { + return x; + } + + __asm__("fabs %d0, %d1" : "=w"(x) : "w"(x)); + return x; +} diff --git a/libm/arch/aarch64/fabsf.c b/libm/arch/aarch64/fabsf.c new file mode 100644 index 00000000..5723764d --- /dev/null +++ b/libm/arch/aarch64/fabsf.c @@ -0,0 +1,30 @@ +/* + * Phoenix-RTOS + * + * libphoenix/libm + * + * Architecture dependent part (arch/aarch64) + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Aleksander Kaminski, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ + +#include + + +float fabsf(float x) +{ + if (isnan(x)) { + return x; + } + if (0 < x) { + return x; + } + + __asm__("fabs %s0, %s1" : "=w"(x) : "w"(x)); + return x; +} diff --git a/libm/arch/aarch64/floord.c b/libm/arch/aarch64/floord.c new file mode 100644 index 00000000..403d4e49 --- /dev/null +++ b/libm/arch/aarch64/floord.c @@ -0,0 +1,27 @@ +/* + * Phoenix-RTOS + * + * libphoenix/libm + * + * Architecture dependent part (arch/aarch64) + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Aleksander Kaminski, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ + +#include + + +double floor(double x) +{ + if (isnan(x) != 0) { + return NAN; + } + + __asm__ volatile("frintm %d0, %d1" : "=w"(x) : "w"(x)); + return x; +} diff --git a/libm/arch/aarch64/floorf.c b/libm/arch/aarch64/floorf.c new file mode 100644 index 00000000..171c1e12 --- /dev/null +++ b/libm/arch/aarch64/floorf.c @@ -0,0 +1,27 @@ +/* + * Phoenix-RTOS + * + * libphoenix/libm + * + * Architecture dependent part (arch/aarch64) + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Aleksander Kaminski, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ + +#include + + +float floorf(float x) +{ + if (isnan(x) != 0) { + return NAN; + } + + __asm__ volatile("frintm %s0, %s1" : "=w"(x) : "w"(x)); + return x; +} diff --git a/libm/arch/aarch64/roundd.c b/libm/arch/aarch64/roundd.c new file mode 100644 index 00000000..ac3ba6e4 --- /dev/null +++ b/libm/arch/aarch64/roundd.c @@ -0,0 +1,27 @@ +/* + * Phoenix-RTOS + * + * libphoenix/libm + * + * Architecture dependent part (arch/aarch64) + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Aleksander Kaminski, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ + +#include + + +double round(double x) +{ + if (isnan(x) != 0) { + return NAN; + } + + __asm__ volatile("frinta %d0, %d1" : "=w"(x) : "w"(x)); + return x; +} diff --git a/libm/arch/aarch64/roundf.c b/libm/arch/aarch64/roundf.c new file mode 100644 index 00000000..dada679d --- /dev/null +++ b/libm/arch/aarch64/roundf.c @@ -0,0 +1,27 @@ +/* + * Phoenix-RTOS + * + * libphoenix/libm + * + * Architecture dependent part (arch/aarch64) + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Aleksander Kaminski, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ + +#include + + +float roundf(float x) +{ + if (isnan(x) != 0) { + return NAN; + } + + __asm__ volatile("frinta %s0, %s1" : "=w"(x) : "w"(x)); + return x; +} diff --git a/libm/arch/aarch64/sqrtd.c b/libm/arch/aarch64/sqrtd.c new file mode 100644 index 00000000..74ea3fc5 --- /dev/null +++ b/libm/arch/aarch64/sqrtd.c @@ -0,0 +1,37 @@ +/* + * Phoenix-RTOS + * + * libphoenix/libm + * + * Architecture dependent part (arch/aarch64) + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Aleksander Kaminski, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ + +#include +#include + + +double sqrt(double x) +{ + if (isnan(x) != 0) { + return NAN; + } + + if (x < 0.0) { + errno = EDOM; + return NAN; + } + + if ((x == 0.0) || (x == INFINITY)) { + return x; + } + + __asm__ volatile("fsqrt %d0, %d1" : "=w"(x) : "w"(x)); + return x; +} diff --git a/libm/arch/aarch64/sqrtf.c b/libm/arch/aarch64/sqrtf.c new file mode 100644 index 00000000..b6fc6f76 --- /dev/null +++ b/libm/arch/aarch64/sqrtf.c @@ -0,0 +1,37 @@ +/* + * Phoenix-RTOS + * + * libphoenix/libm + * + * Architecture dependent part (arch/aarch64) + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Aleksander Kaminski, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ + +#include +#include + + +float sqrtf(float x) +{ + if (isnan(x) != 0) { + return NAN; + } + + if (x < 0.0) { + errno = EDOM; + return NAN; + } + + if ((x == 0.0) || (x == INFINITY)) { + return x; + } + + __asm__ volatile("fsqrt %s0, %s1" : "=w"(x) : "w"(x)); + return x; +} diff --git a/libm/arch/aarch64/truncd.c b/libm/arch/aarch64/truncd.c new file mode 100644 index 00000000..ac310bdf --- /dev/null +++ b/libm/arch/aarch64/truncd.c @@ -0,0 +1,27 @@ +/* + * Phoenix-RTOS + * + * libphoenix/libm + * + * Architecture dependent part (arch/aarch64) + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Aleksander Kaminski, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ + +#include + + +double trunc(double x) +{ + if (isnan(x) != 0) { + return NAN; + } + + __asm__ volatile("frintz %d0, %d1" : "=w"(x) : "w"(x)); + return x; +} diff --git a/libm/arch/aarch64/truncf.c b/libm/arch/aarch64/truncf.c new file mode 100644 index 00000000..b1a90892 --- /dev/null +++ b/libm/arch/aarch64/truncf.c @@ -0,0 +1,27 @@ +/* + * Phoenix-RTOS + * + * libphoenix/libm + * + * Architecture dependent part (arch/aarch64) + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Aleksander Kaminski, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ + +#include + + +float truncf(float x) +{ + if (isnan(x) != 0) { + return NAN; + } + + __asm__ volatile("frintz %s0, %s1" : "=w"(x) : "w"(x)); + return x; +} diff --git a/libm/arch/armv7a/sqrtd.c b/libm/arch/armv7a/sqrtd.c new file mode 100644 index 00000000..31f983df --- /dev/null +++ b/libm/arch/armv7a/sqrtd.c @@ -0,0 +1,36 @@ +/* + * Phoenix-RTOS + * + * libphoenix/libm + * + * Architecture dependent part (arch/armv7a) + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Aleksander Kaminski, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ +#include +#include + + +double sqrt(double x) +{ + if (isnan(x) != 0) { + return NAN; + } + + if (x < 0.0) { + errno = EDOM; + return NAN; + } + + if ((x == 0.0) || (x == INFINITY)) { + return x; + } + + __asm__ volatile("vsqrt.f64 %P0, %P1" : "=w"(x) : "w"(x)); + return x; +} diff --git a/libm/arch/armv7a/sqrtf.c b/libm/arch/armv7a/sqrtf.c new file mode 100644 index 00000000..0b731db6 --- /dev/null +++ b/libm/arch/armv7a/sqrtf.c @@ -0,0 +1,36 @@ +/* + * Phoenix-RTOS + * + * libphoenix/libm + * + * Architecture dependent part (arch/armv7a) + * + * Copyright 2017, 2018, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Aleksander Kaminski, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ +#include +#include + + +float sqrtf(float x) +{ + if (isnan(x) != 0) { + return NAN; + } + + if (x < 0.0) { + errno = EDOM; + return NAN; + } + + if ((x == 0.0) || (x == INFINITY)) { + return x; + } + + __asm__ volatile("vsqrt.f32 %0, %1" : "=t"(x) : "t"(x)); + return x; +} diff --git a/libm/arch/armv7m/.gitkeep b/libm/arch/armv7m/.gitkeep new file mode 100644 index 00000000..e69de29b diff --git a/libm/arch/armv7r/sqrtd.c b/libm/arch/armv7r/sqrtd.c new file mode 100644 index 00000000..3b0ee0d4 --- /dev/null +++ b/libm/arch/armv7r/sqrtd.c @@ -0,0 +1,36 @@ +/* + * Phoenix-RTOS + * + * libphoenix/libm + * + * Architecture dependent part (arch/armv7r) + * + * Copyright 2017, 2024, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Lukasz Leczkowski, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ +#include +#include + + +double sqrt(double x) +{ + if (isnan(x) != 0) { + return NAN; + } + + if (x < 0.0) { + errno = EDOM; + return NAN; + } + + if ((x == 0.0) || (x == INFINITY)) { + return x; + } + + __asm__ volatile("vsqrt.f64 %P0, %P1" : "=w"(x) : "w"(x)); + return x; +} diff --git a/libm/arch/armv7r/sqrtf.c b/libm/arch/armv7r/sqrtf.c new file mode 100644 index 00000000..261b8560 --- /dev/null +++ b/libm/arch/armv7r/sqrtf.c @@ -0,0 +1,36 @@ +/* + * Phoenix-RTOS + * + * libphoenix/libm + * + * Architecture dependent part (arch/armv7r) + * + * Copyright 2017, 2024, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Lukasz Leczkowski, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ +#include +#include + + +float sqrtf(float x) +{ + if (isnan(x) != 0) { + return NAN; + } + + if (x < 0.0) { + errno = EDOM; + return NAN; + } + + if ((x == 0.0) || (x == INFINITY)) { + return x; + } + + __asm__ volatile("vsqrt.f32 %0, %1" : "=t"(x) : "t"(x)); + return x; +} diff --git a/libm/arch/armv8m/.gitkeep b/libm/arch/armv8m/.gitkeep new file mode 100644 index 00000000..e69de29b diff --git a/libm/arch/armv8r/sqrtd.c b/libm/arch/armv8r/sqrtd.c new file mode 100644 index 00000000..622cfc22 --- /dev/null +++ b/libm/arch/armv8r/sqrtd.c @@ -0,0 +1,37 @@ +/* + * Phoenix-RTOS + * + * libphoenix/libm + * + * Architecture dependent part (arch/armv8r) + * + * Copyright 2017, 2024, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Lukasz Leczkowski, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ + +#include +#include + + +double sqrt(double x) +{ + if (isnan(x) != 0) { + return NAN; + } + + if (x < 0.0) { + errno = EDOM; + return NAN; + } + + if ((x == 0.0) || (x == INFINITY)) { + return x; + } + + __asm__ volatile("vsqrt.f64 %P0, %P1" : "=w"(x) : "w"(x)); + return x; +} diff --git a/libm/arch/armv8r/sqrtf.c b/libm/arch/armv8r/sqrtf.c new file mode 100644 index 00000000..62c0a7d7 --- /dev/null +++ b/libm/arch/armv8r/sqrtf.c @@ -0,0 +1,37 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * Architecture dependent part (arch/armv8r) + * + * Copyright 2017, 2024, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Lukasz Leczkowski, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ + +#include +#include + + +float sqrtf(float x) +{ + if (isnan(x) != 0) { + return NAN; + } + + if (x < 0.0) { + errno = EDOM; + return NAN; + } + + if ((x == 0.0) || (x == INFINITY)) { + return x; + } + + __asm__ volatile("vsqrt.f32 %0, %1" : "=t"(x) : "t"(x)); + return x; +} diff --git a/libm/arch/ia32/sqrtd.c b/libm/arch/ia32/sqrtd.c new file mode 100644 index 00000000..1e3923cc --- /dev/null +++ b/libm/arch/ia32/sqrtd.c @@ -0,0 +1,54 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * Architecture dependent part + * + * Copyright 2017, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ + +#include +#include + + +double sqrt(double x) +{ + if (isnan(x) != 0) { + return NAN; + } + + if (x < 0.0) { + errno = EDOM; + return NAN; + } + + if ((x == 0.0) || (x == INFINITY)) { + return x; + } + + unsigned short newcw, savecw; + unsigned short register reg; + + /* clang-format off */ + __asm__ volatile ( + "fstcw %w[savecw]\n\t" /* save control word */ + "mov %w[savecw], %w[reg]\n\t" /* reg <- (memptr) */ + "and $0xf0ff, %w[reg]\n\t" /* inherit exception mask and infinity */ + "or $0x200, %w[reg]\n\t" /* set double precision, round to nearest-even */ + "mov %w[reg], %w[newcw]\n\t" /* (memptr) <- reg */ + "fldcw %w[newcw]\n\t" /* and apply (memptr) new control word */ + "fsqrt\n\t" /* calculate double precision sqrt */ + "fldcw %w[savecw]" /* restore control word */ + : "+t" (x), [reg] "=&r" (reg) + : [newcw] "m" (newcw), [savecw] "m" (savecw) + : "cc", "memory"); + /* clang-format on */ + + return x; +} diff --git a/libm/arch/ia32/sqrtf.c b/libm/arch/ia32/sqrtf.c new file mode 100644 index 00000000..9b046adf --- /dev/null +++ b/libm/arch/ia32/sqrtf.c @@ -0,0 +1,36 @@ +/* + * Phoenix-RTOS + * + * libphoenix + * + * Architecture dependent part + * + * Copyright 2017, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ + +#include +#include + + +float sqrtf(float x) +{ + if (isnan(x) != 0) { + return NAN; + } + + if (x < 0.0) { + errno = EDOM; + return NAN; + } + + if ((x == 0.0) || (x == INFINITY)) { + return x; + } + + return (float)sqrt((float)x); +} diff --git a/libm/arch/riscv64/sqrtd.c b/libm/arch/riscv64/sqrtd.c new file mode 100644 index 00000000..8d8096b6 --- /dev/null +++ b/libm/arch/riscv64/sqrtd.c @@ -0,0 +1,39 @@ +/* + * Phoenix-RTOS + * + * libphoenix/libm + * + * Architecture dependent part (RISCV64) + * + * Copyright 2018, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ + +#include +#include + +double sqrt(double x) +{ + if (isnan(x) != 0) { + return NAN; + } + + if (x < 0.0) { + errno = EDOM; + return NAN; + } + + if ((x == 0.0) || (x == INFINITY)) { + return x; + } + + /* clang-format off */ + __asm__ volatile ("fsqrt.d %0, %1" : "=f"(x) : "f"(x)); + /* clang-format on */ + + return x; +} diff --git a/libm/arch/riscv64/sqrtf.c b/libm/arch/riscv64/sqrtf.c new file mode 100644 index 00000000..332b55b3 --- /dev/null +++ b/libm/arch/riscv64/sqrtf.c @@ -0,0 +1,40 @@ +/* + * Phoenix-RTOS + * + * libphoenix/libm + * + * Architecture dependent part (RISCV64) + * + * Copyright 2018, 2026 Phoenix Systems + * Author: Pawel Pisarczyk, Michal Lach + * + * This file is part of Phoenix-RTOS. + * + * %LICENSE% + */ + +#include +#include + + +float sqrtf(float x) +{ + if (isnan(x) != 0) { + return NAN; + } + + if (x < 0.0) { + errno = EDOM; + return NAN; + } + + if ((x == 0.0) || (x == INFINITY)) { + return x; + } + + /* clang-format off */ + __asm__ volatile ("fsqrt.s %0, %1" : "=f"(x) : "f"(x)); + /* clang-format on */ + + return x; +} diff --git a/libm/arch/sparcv8leon/.gitkeep b/libm/arch/sparcv8leon/.gitkeep new file mode 100644 index 00000000..e69de29b