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