michaelh/zephyr
1
0
Fork 0
Code Releases Activity

libc: minimal: math: Removed undefined behavior in sqrt routines

Browse source 
The previous code used casts of address of int to float pointer.
This is undefined behavior in C and may cause a compiler to
optimize code incorrectly.

I have replaced it with a union of int and float types,
eliminating the cast problem.

Signed-off-by: Lars-Ove Karlsson <lars-ove.karlsson@iar.com>
This commit is contained in:
Lars-Ove Karlsson 2024-09-26 14:28:04 +02:00 • committed by Fabio Baltieri
commit 0a8e1ad972
2 changed files with 32 additions and 24 deletions
Download patch file Download diff file Expand all files Collapse all files

28
lib/libc/minimal/source/math/sqrt.c
View file

@ -15,15 +15,19 @@
#define MAX_D_ERROR_COUNT 5 /* when result almost converges, stop */
#define EXP_MASK64 GENMASK64(62, 52)
typedef union {
double d;
int64_t i;
} int64double_t;
double sqrt(double square)
{
int i;
int64_t exponent;
double root;
double last;
int64_t *p_square = (int64_t *)&square;
int64_t *p_root = (int64_t *)&root;
int64_t *p_last = (int64_t *)&last;
int64_t exponent;
int64double_t root;
int64double_t last;
int64double_t p_square;
p_square.d = square;
if (square == 0.0) {
return square;
@ -36,21 +40,21 @@ double sqrt(double square)
* we can do this by dividing the exponent part of the float by 2
* this assumes IEEE-754 format doubles
*/
exponent = ((*p_square & EXP_MASK64)>>52)-1023;
if (exponent == 0x7FF-1023) {
exponent = ((p_square.i & EXP_MASK64) >> 52) - 1023;
if (exponent == 0x7FF - 1023) {
/* the number is a NAN or inf, return NaN or inf */
return square + square;
}
exponent /= 2;
*p_root = (*p_square & ~EXP_MASK64) | (exponent+1023)<<52;
root.i = (p_square.i & ~EXP_MASK64) | (exponent + 1023) << 52;
for (i = 0; i < MAX_D_ITTERATIONS; i++) {
last = root;
root = (root + square / root) * 0.5;
root.d = (root.d + square / root.d) * 0.5;
/* if (llabs(*p_root-*p_last)<MAX_D_ERROR_COUNT) */
if ((*p_root ^ *p_last) < MAX_D_ERROR_COUNT) {
if ((root.i ^ last.i) < MAX_D_ERROR_COUNT) {
break;
}
}
return root;
return root.d;
}

28
lib/libc/minimal/source/math/sqrtf.c
View file

@ -15,15 +15,19 @@
#define MAX_F_ERROR_COUNT 3 /* when result almost converges, stop */
#define EXP_MASK32 GENMASK(30, 23)
typedef union {
float f;
int32_t i;
} intfloat_t;
float sqrtf(float square)
{
int i;
float root;
float last;
int32_t exponent;
int32_t *p_square = (int32_t *)&square;
int32_t *p_root = (int32_t *)&root;
int32_t *p_last = (int32_t *)&last;
int32_t exponent;
intfloat_t root;
intfloat_t last;
intfloat_t p_square;
p_square.f = square;
if (square == 0.0f) {
return square;
@ -36,21 +40,21 @@ float sqrtf(float square)
* we can do this by dividing the exponent part of the float by 2
* this assumes IEEE-754 format doubles
*/
exponent = ((*p_square & EXP_MASK32)>>23)-127;
if (exponent == 0xFF-127) {
exponent = ((p_square.i & EXP_MASK32) >> 23) - 127;
if (exponent == 0xFF - 127) {
/* the number is a NAN or inf, return NaN or inf */
return square + square;
}
exponent /= 2;
*p_root = (*p_square & ~EXP_MASK32) | (exponent+127) << 23;
root.i = (p_square.i & ~EXP_MASK32) | (exponent + 127) << 23;
for (i = 0; i < MAX_F_ITTERATIONS; i++) {
last = root;
root = (root + square / root) * 0.5f;
root.f = (root.f + square / root.f) * 0.5f;
/* if (labs(*p_root - *p_last) < MAX_F_ERROR_COUNT) */
if ((*p_root ^ *p_last) < MAX_F_ERROR_COUNT) {
if ((root.i ^ last.i) < MAX_F_ERROR_COUNT) {
break;
}
}
return root;
return root.f;
}