1 /*
2  * Copyright (c) 2023 Lawrence King
3  *
4  * SPDX-License-Identifier: Apache-2.0
5  *
6  */
7 
8 #include <math.h>
9 #include <zephyr/kernel.h>
10 #include <zephyr/ztest.h>
11 
12 #define local_abs(x) (((x) < 0) ? -(x) : (x))
13 
14 #ifndef NAN
15 #define NAN	(__builtin_nansf(""))
16 #endif
17 
18 #ifndef INFINITY
19 #define INFINITY	(__builtin_inff())
20 #endif
21 
22 static float test_floats[] = {
23 	1.0f, 2.0f, 3.0f, 4.0f,
24 	5.0f, 6.0f, 7.0f, 8.0f, 9.0f,	/* numbers across the decade */
25 	3.14159265359f, 2.718281828f,	/* irrational numbers pi and e */
26 	123.4f,	0.025f,	0.10f,	1.875f	/* numbers with infinite */
27 					/* repeating binary representation */
28 	};
29 #define	NUM_TEST_FLOATS	(sizeof(test_floats)/sizeof(float))
30 
31 static double test_doubles[] = {
32 	1.0, 2.0, 3.0, 4.0,
33 	5.0, 6.0, 7.0, 8.0, 9.0,	/* numbers across the decade */
34 	3.14159265359, 2.718281828,	/* irrational numbers pi and e */
35 	123.4,	0.025,	0.10,	1.875	/* numbers with infinite */
36 					/* repeating binary representationa */
37 };
38 #define	NUM_TEST_DOUBLES	(sizeof(test_floats)/sizeof(float))
39 
40 #ifndef	isinf
isinf(double x)41 static int isinf(double x)
42 {
43 	union { uint64_t u; double d; } ieee754;
44 	ieee754.d = x;
45 	ieee754.u &= ~0x8000000000000000;	/* ignore the sign */
46 	return ((ieee754.u >> 52) == 0x7FF) &&
47 		((ieee754.u & 0x000fffffffffffff) == 0);
48 }
49 #endif
50 
51 #ifndef	isnan
isnan(double x)52 static int isnan(double x)
53 {
54 	union { uint64_t u; double d; } ieee754;
55 	ieee754.d = x;
56 	ieee754.u &= ~0x8000000000000000;	/* ignore the sign */
57 	return ((ieee754.u >> 52) == 0x7FF) &&
58 		((ieee754.u & 0x000fffffffffffff) != 0);
59 }
60 #endif
61 
62 #ifndef	isinff
isinff(float x)63 static int isinff(float x)
64 {
65 	union { uint32_t u; float f; } ieee754;
66 	ieee754.f = x;
67 	ieee754.u &= ~0x80000000;		/* ignore the sign */
68 	return ((ieee754.u >> 23) == 0xFF) &&
69 		((ieee754.u & 0x7FFFFF) == 0);
70 }
71 #endif
72 
73 #ifndef	isnanf
isnanf(float x)74 static int isnanf(float x)
75 {
76 	union { uint32_t u; float f; } ieee754;
77 	ieee754.f = x;
78 	ieee754.u &= ~0x80000000;		/* ignore the sign */
79 	return ((ieee754.u >> 23) == 0xFF) &&
80 		((ieee754.u & 0x7FFFFF) != 0);
81 }
82 #endif
83 
84 /* small errors are expected, computed as percentage error */
85 #define MAX_FLOAT_ERROR_PERCENT		(3.5e-5f)
86 #define MAX_DOUBLE_ERROR_PERCENT	(4.5e-14)
87 
ZTEST(libc_common,test_sqrtf)88 ZTEST(libc_common, test_sqrtf)
89 {
90 	int i;
91 	float exponent, resf, square, root_squared, error;
92 	uint32_t max_error;
93 	int32_t ierror;
94 	int32_t *p_square = (int32_t *)&square;
95 	int32_t *p_root_squared = (int32_t *)&root_squared;
96 
97 	max_error = 0;
98 
99 	/* test the special cases of 0.0, NAN, -NAN, INFINITY, -INFINITY,  and -10.0 */
100 	zassert_true(sqrtf(0.0f) == 0.0f, "sqrtf(0.0)");
101 	zassert_true(isnanf(sqrtf(NAN)), "sqrt(nan)");
102 #ifdef issignallingf
103 	zassert_true(issignallingf(sqrtf(NAN)), "ssignalingf(sqrtf(nan))");
104 /*	printf("issignallingf();\n"); */
105 #endif
106 	zassert_true(isnanf(sqrtf(-NAN)), "isnanf(sqrtf(-nan))");
107 	zassert_true(isinff(sqrtf(INFINITY)), "isinff(sqrt(inf))");
108 	zassert_true(isnanf(sqrtf(-INFINITY)), "isnanf(sqrt(-inf))");
109 	zassert_true(isnanf(sqrtf(-10.0f)), "isnanf(sqrt(-10.0))");
110 
111 	for (exponent = 1.0e-10f; exponent < 1.0e10f; exponent *= 10.0f) {
112 		for (i = 0; i < NUM_TEST_FLOATS; i++) {
113 			square = test_floats[i] * exponent;
114 			resf = sqrtf(square);
115 			root_squared = resf * resf;
116 			zassert_true((resf > 0.0f) && (resf < INFINITY),
117 					"sqrtf out of range");
118 			if ((resf > 0.0f) && (resf < INFINITY)) {
119 				error = (square - root_squared) /
120 					square * 100;
121 				if (error < 0.0f) {
122 					error = -error;
123 				}
124 				/* square and root_squared should be almost identical
125 				 * except the last few bits, the EXOR will only set
126 				 * the bits that are different
127 				 */
128 				ierror = (*p_square - *p_root_squared);
129 				ierror = local_abs(ierror);
130 				if (ierror > max_error) {
131 					max_error = ierror;
132 				}
133 			} else {
134 				/* negative, +NaN, -NaN, inf or -inf */
135 				error = 0.0f;
136 			}
137 			zassert_true(error < MAX_FLOAT_ERROR_PERCENT,
138 					"max sqrtf error exceeded");
139 		}
140 	}
141 	zassert_true(max_error < 0x03, "huge errors in sqrt implementation");
142 	/* print the max error */
143 	TC_PRINT("test_sqrtf max error %d counts\n", max_error);
144 }
145 
ZTEST(libc_common,test_sqrt)146 ZTEST(libc_common, test_sqrt)
147 {
148 	if (sizeof(double) != 8) {
149 		ztest_test_skip();
150 	}
151 
152 	int i;
153 	double resd, error, square, root_squared, exponent;
154 	uint64_t max_error;
155 	int64_t ierror;
156 	int64_t *p_square = (int64_t *)&square;
157 	int64_t *p_root_squared = (int64_t *)&root_squared;
158 
159 	max_error = 0;
160 
161 	/* test the special cases of 0.0, NAN, -NAN, INFINITY, -INFINITY,  and -10.0 */
162 	zassert_true(sqrt(0.0) == 0.0, "sqrt(0.0)");
163 	zassert_true(isnan(sqrt((double)NAN)), "sqrt(nan)");
164 #ifdef issignalling
165 	zassert_true(issignalling(sqrt((double)NAN)), "ssignaling(sqrt(nan))");
166 /*	printf("issignalling();\n"); */
167 #endif
168 	zassert_true(isnan(sqrt((double)-NAN)), "isnan(sqrt(-nan))");
169 	zassert_true(isinf(sqrt((double)INFINITY)), "isinf(sqrt(inf))");
170 	zassert_true(isnan(sqrt((double)-INFINITY)), "isnan(sqrt(-inf))");
171 	zassert_true(isnan(sqrt(-10.0)), "isnan(sqrt(-10.0))");
172 
173 	for (exponent = 1.0e-10; exponent < 1.0e10; exponent *= 10.0) {
174 		for (i = 0; i < NUM_TEST_DOUBLES; i++) {
175 			square = test_doubles[i] * exponent;
176 			resd = sqrt(square);
177 			root_squared = resd * resd;
178 			zassert_true((resd > 0.0) && (resd < (double)INFINITY),
179 					"sqrt out of range");
180 			if ((resd > 0.0) && (resd < (double)INFINITY)) {
181 				error = (square - root_squared) /
182 					square * 100;
183 				if (error < 0.0) {
184 					error = -error;
185 				}
186 				/* square and root_squared should be almost identical
187 				 * except the last few bits, the EXOR will only set
188 				 * the bits that are different
189 				 */
190 				ierror = (*p_square - *p_root_squared);
191 				ierror = local_abs(ierror);
192 				if (ierror > max_error) {
193 					max_error = ierror;
194 				}
195 			} else {
196 				/* negative, +NaN, -NaN, inf or -inf */
197 				error = 0.0;
198 			}
199 			zassert_true(error < MAX_DOUBLE_ERROR_PERCENT,
200 					"max sqrt error exceeded");
201 		}
202 	}
203 	zassert_true(max_error < 0x04, "huge errors in sqrt implementation");
204 	/* print the max error */
205 	TC_PRINT("test_sqrt max error %d counts\n", (uint32_t)max_error);
206 }
207