1 #include <iostream>
2 #include <iomanip>
3 #include <limits>
4 #include <type_traits>
5 #include "vnl/vnl_math.h"
6 #include "vnl/vnl_complex.h" // for vnl_math::abs(std::complex)
7 #include "testlib/testlib_test.h"
8
9 static constexpr double vnl_math_test_20_epsilon = 2 * 10 * std::numeric_limits<double>::epsilon();
10 static constexpr double vnl_math_test_2_epsilon = 2 * std::numeric_limits<double>::epsilon();
11
12 // Utility function for printing hex representations
13 template <typename T>
14 std::string
print_hex(const T p)15 print_hex(const T p)
16 {
17 std::stringstream str;
18 str << std::hex << std::setfill('0') << std::setw(2);
19 for (int i = 0; i < static_cast<int>(16 - sizeof(p)); ++i)
20 {
21 str << "..";
22 }
23 for (int i = (sizeof(p) - 1); i >= 0; --i)
24 {
25 str << std::setfill('0') << std::setw(2);
26 const auto curr_value = static_cast<short>((reinterpret_cast<unsigned char const *>(&p))[i]);
27 str << curr_value;
28 }
29 str << std::dec;
30 return str.str();
31 }
32
33 static void
check_pointer(const void *)34 check_pointer(const void *)
35 {}
36
37 static void
test_static_const_definition()38 test_static_const_definition()
39 {
40 check_pointer(&vnl_math::e);
41 check_pointer(&vnl_math::euler);
42 check_pointer(&vnl_math::log2e);
43 check_pointer(&vnl_math::log10e);
44 check_pointer(&vnl_math::ln2);
45 check_pointer(&vnl_math::ln10);
46 check_pointer(&vnl_math::pi);
47 check_pointer(&vnl_math::twopi);
48 check_pointer(&vnl_math::pi_over_2);
49 check_pointer(&vnl_math::pi_over_4);
50 check_pointer(&vnl_math::pi_over_180);
51 check_pointer(&vnl_math::one_over_pi);
52 check_pointer(&vnl_math::two_over_pi);
53 check_pointer(&vnl_math::sqrt2pi);
54 check_pointer(&vnl_math::one_over_sqrt2pi);
55 check_pointer(&vnl_math::two_over_sqrtpi);
56 check_pointer(&vnl_math::deg_per_rad);
57 check_pointer(&vnl_math::sqrt2);
58 check_pointer(&vnl_math::sqrt1_2);
59 check_pointer(&vnl_math::sqrt1_3);
60 check_pointer(&vnl_math::eps);
61 check_pointer(&vnl_math::sqrteps);
62 }
63
64 // Test that the vnl_math constants don't have weird values
65 static void
test_math_constants()66 test_math_constants()
67 {
68 #define TEST_CONSTANT(a, v) TEST_NEAR("value: ", vnl_math::a, v, 0);
69 TEST_NEAR("log of e is 1", log(vnl_math::e), 1.0, 1e-15);
70 TEST_CONSTANT(e, 2.7182818284590452353602874713526624977572470936999);
71 TEST_NEAR("log2e * ln2 = 1", vnl_math::log2e * vnl_math::ln2, 1.0, 1e-15);
72 TEST_CONSTANT(log2e, 1.4426950408889634073599246810018921374266459541529);
73 TEST_CONSTANT(ln2, 0.69314718055994530941723212145817656807550013436025);
74 TEST_NEAR("log10e * ln10 = 1", vnl_math::log10e * vnl_math::ln10, 1.0, 1e-15);
75 TEST_CONSTANT(log10e, 0.43429448190325182765112891891660508229439700580366);
76 TEST_CONSTANT(ln10, 2.3025850929940456840179914546843642076011014886287);
77 TEST_NEAR("cos(pi) = -1", cos(vnl_math::pi), -1.0, 1e-15);
78 TEST_CONSTANT(pi, 3.1415926535897932384626433832795028841971693993751);
79 TEST_NEAR("twopi = 2*pi", vnl_math::twopi, 2.0 * vnl_math::pi, 1e-15);
80 TEST_CONSTANT(twopi, 6.2831853071795864769252867665590057683943387987502);
81 TEST_NEAR("pi_over_2 = pi/2", vnl_math::pi_over_2, 0.5 * vnl_math::pi, 1e-15);
82 TEST_CONSTANT(pi_over_2, 1.5707963267948966192313216916397514420985846996875);
83 TEST_NEAR("pi_over_4 = pi/4", vnl_math::pi_over_4, 0.25 * vnl_math::pi, 1e-15);
84 TEST_CONSTANT(pi_over_4, 0.78539816339744830961566084581987572104929234984377);
85 TEST_NEAR("pi_over_180=pi/180", vnl_math::pi_over_180, vnl_math::pi / 180.0, 1e-15);
86 TEST_CONSTANT(pi_over_180, 0.017453292519943295769236907684886127134428718885417);
87 TEST_NEAR("pi*one_over_pi=1", vnl_math::pi * vnl_math::one_over_pi, 1.0, 1e-15);
88 TEST_CONSTANT(one_over_pi, 0.31830988618379067153776752674502872406891929148091);
89 TEST_NEAR("pi*two_over_pi=2", vnl_math::pi * vnl_math::two_over_pi, 2.0, 1e-15);
90 TEST_CONSTANT(two_over_pi, 0.63661977236758134307553505349005744813783858296182);
91 TEST_NEAR("deg_per_rad=180/pi", vnl_math::deg_per_rad, 180.0 / vnl_math::pi, 1e-15);
92 TEST_CONSTANT(deg_per_rad, 57.295779513082320876798154814105170332405472466564);
93 TEST_NEAR("sqrt2pi^2", vnl_math::sqrt2pi * vnl_math::sqrt2pi, vnl_math::twopi, 1e-15);
94 TEST_CONSTANT(sqrt2pi, 2.5066282746310005024157652848110452530069867406099);
95 TEST_NEAR("two_over_sqrtpi", vnl_math::two_over_sqrtpi, 2.0 / sqrt(vnl_math::pi), 1e-15);
96 TEST_CONSTANT(two_over_sqrtpi, 1.1283791670955125738961589031215451716881012586579);
97 TEST_NEAR("one_over_sqrt2pi", vnl_math::one_over_sqrt2pi, 1.0 / sqrt(vnl_math::twopi), 1e-15);
98 TEST_CONSTANT(one_over_sqrt2pi, 0.39894228040143267793994605993438186847585863116493);
99 TEST_NEAR("sqrt2*sqrt2=2", vnl_math::sqrt2 * vnl_math::sqrt2, 2.0, 1e-15);
100 TEST_CONSTANT(sqrt2, 1.4142135623730950488016887242096980785696718753769);
101 TEST_NEAR("sqrt1_2*sqrt2=1", vnl_math::sqrt1_2 * vnl_math::sqrt2, 1.0, 1e-15);
102 TEST_CONSTANT(sqrt1_2, 0.70710678118654752440084436210484903928483593768847);
103 TEST_NEAR("sqrt1_3^2=1/3", vnl_math::sqrt1_3 * vnl_math::sqrt1_3, 1.0 / 3.0, 1e-15);
104 TEST_CONSTANT(sqrt1_3, 0.57735026918962576450914878050195745564760175127012);
105 TEST_CONSTANT(euler, 0.57721566490153286060651209008240243104215933593992);
106 #undef TEST_CONSTANT
107 }
108
109 static void
test_math()110 test_math()
111 {
112 // Call it to avoid compiler warnings
113 test_static_const_definition();
114 test_math_constants();
115
116 int n = -11;
117 float f = -7.5f;
118 double d = -vnl_math::pi;
119 std::complex<double> i(0, 1);
120 std::complex<double> z(-1, 2);
121 std::complex<double> e_ipi = std::exp(d * i);
122
123 std::cout << "n = " << n << '\n'
124 << "f = " << f << '\n'
125 << "d = " << d << '\n'
126 << "i = " << i << '\n'
127 << "z = " << z << '\n'
128 << "exp(d*i) = " << e_ipi << '\n'
129 << '\n'
130
131 << "abs(n) = " << vnl_math::abs(n) << '\n'
132 << "abs(f) = " << vnl_math::abs(f) << '\n'
133 << "abs(d) = " << vnl_math::abs(d) << '\n'
134 << "abs(i) = " << vnl_math::abs(i) << '\n'
135 << "abs(z) = " << vnl_math::abs(z) << '\n'
136 << "norm(z) = " << vnl_math::squared_magnitude(z) << '\n'
137 << std::endl;
138
139 TEST("abs(n) == 11", vnl_math::abs(n), 11);
140 TEST("abs(f) == 7.5f", vnl_math::abs(f), 7.5f);
141 TEST("abs(d) == pi", vnl_math::abs(d), vnl_math::pi);
142 TEST("abs(i) == 1", vnl_math::abs(i), 1.0);
143 TEST_NEAR("abs(-1+2i)~=sqrt(5)", vnl_math::abs(z), std::sqrt(5.0), 1e-12);
144 TEST_NEAR("norm(-1+2i) ~= 5", vnl_math::squared_magnitude(z), 5, 1e-12);
145 TEST_NEAR("exp(d*i) ~= -1", vnl_math::abs(e_ipi + 1.0), 0, 1e-12);
146 std::cout << std::endl;
147
148 TEST("rnd(-8.4999) == -8 ", vnl_math::rnd(-8.4999), -8);
149 TEST("rnd(-8.4999f) == -8 ", vnl_math::rnd(-8.4999f), -8);
150 TEST("rnd(-8.50) == -8/9", vnl_math::rnd(-8.50) / 2, -4);
151 TEST("rnd(-8.50f) == -8/9", vnl_math::rnd(-8.50f) / 2, -4);
152 TEST("rnd(-8.5001) == -9 ", vnl_math::rnd(-8.5001), -9);
153 TEST("rnd(-8.5001f) == -9 ", vnl_math::rnd(-8.5001f), -9);
154 TEST("rnd(8.4999) == 8 ", vnl_math::rnd(8.4999), 8);
155 TEST("rnd(8.4999f) == 8 ", vnl_math::rnd(8.4999f), 8);
156 TEST("rnd(8.50) == 8/9", vnl_math::rnd(8.50) / 2, 4);
157 TEST("rnd(8.50f) == 8/9", vnl_math::rnd(8.50f) / 2, 4);
158 TEST("rnd(8.5001) == 9 ", vnl_math::rnd(8.5001), 9);
159 TEST("rnd(8.5001f) == 9 ", vnl_math::rnd(8.5001f), 9);
160
161 TEST("rnd(-9.4999) == -9 ", vnl_math::rnd(-9.4999), -9);
162 TEST("rnd(-9.4999f) == -9 ", vnl_math::rnd(-9.4999f), -9);
163 TEST("rnd(-9.50) == -9/10", (vnl_math::rnd(-9.50) + 1) / 2, -4);
164 TEST("rnd(-9.50f) == -9/10", (vnl_math::rnd(-9.50f) + 1) / 2, -4);
165 TEST("rnd(-9.5001) == -10 ", vnl_math::rnd(-9.5001), -10);
166 TEST("rnd(-9.5001f) == -10 ", vnl_math::rnd(-9.5001f), -10);
167 TEST("rnd(9.4999) == 9 ", vnl_math::rnd(9.4999), 9);
168 TEST("rnd(9.4999f) == 9 ", vnl_math::rnd(9.4999f), 9);
169 TEST("rnd(9.50) == 9/10", (vnl_math::rnd(9.50) - 1) / 2, 4);
170 TEST("rnd(9.50f) == 9/10", (vnl_math::rnd(9.50f) - 1) / 2, 4);
171 TEST("rnd(9.5001) == 10 ", vnl_math::rnd(9.5001), 10);
172 TEST("rnd(9.5001f) == 10 ", vnl_math::rnd(9.5001f), 10);
173
174 TEST("rnd_halfinttoeven(-8.4999) == -8", vnl_math::rnd_halfinttoeven(-8.4999), -8);
175 TEST("rnd_halfinttoeven(-8.4999f) == -8", vnl_math::rnd_halfinttoeven(-8.4999f), -8);
176 TEST("rnd_halfinttoeven(-8.50) == -8", vnl_math::rnd_halfinttoeven(-8.50), -8);
177 TEST("rnd_halfinttoeven(-8.50f) == -8", vnl_math::rnd_halfinttoeven(-8.50f), -8);
178 TEST("rnd_halfinttoeven(-8.5001) == -9", vnl_math::rnd_halfinttoeven(-8.5001), -9);
179 TEST("rnd_halfinttoeven(-8.5001f) == -9", vnl_math::rnd_halfinttoeven(-8.5001f), -9);
180 TEST("rnd_halfinttoeven(8.4999) == 8", vnl_math::rnd_halfinttoeven(8.4999), 8);
181 TEST("rnd_halfinttoeven(8.4999f) == 8", vnl_math::rnd_halfinttoeven(8.4999f), 8);
182 TEST("rnd_halfinttoeven(8.50) == 8", vnl_math::rnd_halfinttoeven(8.50), 8);
183 TEST("rnd_halfinttoeven(8.50f) == 8", vnl_math::rnd_halfinttoeven(8.50f), 8);
184 TEST("rnd_halfinttoeven(8.5001) == 9", vnl_math::rnd_halfinttoeven(8.5001), 9);
185 TEST("rnd_halfinttoeven(8.5001f) == 9", vnl_math::rnd_halfinttoeven(8.5001f), 9);
186
187 TEST("rnd_halfinttoeven(-9.4999) == -9 ", vnl_math::rnd_halfinttoeven(-9.4999), -9);
188 TEST("rnd_halfinttoeven(-9.4999f) == -9 ", vnl_math::rnd_halfinttoeven(-9.4999f), -9);
189 TEST("rnd_halfinttoeven(-9.50) == -10", vnl_math::rnd_halfinttoeven(-9.50), -10);
190 TEST("rnd_halfinttoeven(-9.50f) == -10", vnl_math::rnd_halfinttoeven(-9.50f), -10);
191 TEST("rnd_halfinttoeven(-9.5001) == -10", vnl_math::rnd_halfinttoeven(-9.5001), -10);
192 TEST("rnd_halfinttoeven(-9.5001f) == -10", vnl_math::rnd_halfinttoeven(-9.5001f), -10);
193 TEST("rnd_halfinttoeven(9.4999) == 9 ", vnl_math::rnd_halfinttoeven(9.4999), 9);
194 TEST("rnd_halfinttoeven(9.4999f) == 9 ", vnl_math::rnd_halfinttoeven(9.4999f), 9);
195 TEST("rnd_halfinttoeven(9.50) == 10", vnl_math::rnd_halfinttoeven(9.50), 10);
196 TEST("rnd_halfinttoeven(9.50f) == 10", vnl_math::rnd_halfinttoeven(9.50f), 10);
197 TEST("rnd_halfinttoeven(9.5001) == 10", vnl_math::rnd_halfinttoeven(9.5001), 10);
198 TEST("rnd_halfinttoeven(9.5001f) == 10", vnl_math::rnd_halfinttoeven(9.5001f), 10);
199
200 TEST("rnd_halfintup(-8.4999) == -8", vnl_math::rnd_halfintup(-8.4999), -8);
201 TEST("rnd_halfintup(-8.4999f) == -8", vnl_math::rnd_halfintup(-8.4999f), -8);
202 TEST("rnd_halfintup(-8.50) == -8", vnl_math::rnd_halfintup(-8.50), -8);
203 TEST("rnd_halfintup(-8.50f) == -8", vnl_math::rnd_halfintup(-8.50f), -8);
204 TEST("rnd_halfintup(-8.5001) == -9", vnl_math::rnd_halfintup(-8.5001), -9);
205 TEST("rnd_halfintup(-8.5001f) == -9", vnl_math::rnd_halfintup(-8.5001f), -9);
206 TEST("rnd_halfintup(8.4999) == 8", vnl_math::rnd_halfintup(8.4999), 8);
207 TEST("rnd_halfintup(8.4999f) == 8", vnl_math::rnd_halfintup(8.4999f), 8);
208 TEST("rnd_halfintup(8.50) == 9", vnl_math::rnd_halfintup(8.50), 9);
209 TEST("rnd_halfintup(8.50f) == 9", vnl_math::rnd_halfintup(8.50f), 9);
210 TEST("rnd_halfintup(8.5001) == 9", vnl_math::rnd_halfintup(8.5001), 9);
211 TEST("rnd_halfintup(8.5001f) == 9", vnl_math::rnd_halfintup(8.5001f), 9);
212
213 TEST("rnd_halfintup(-9.4999) == -9 ", vnl_math::rnd_halfintup(-9.4999), -9);
214 TEST("rnd_halfintup(-9.4999f) == -9 ", vnl_math::rnd_halfintup(-9.4999f), -9);
215 TEST("rnd_halfintup(-9.50) == -9 ", vnl_math::rnd_halfintup(-9.50), -9);
216 TEST("rnd_halfintup(-9.50f) == -9 ", vnl_math::rnd_halfintup(-9.50f), -9);
217 TEST("rnd_halfintup(-9.5001) == -10", vnl_math::rnd_halfintup(-9.5001), -10);
218 TEST("rnd_halfintup(-9.5001f) == -10", vnl_math::rnd_halfintup(-9.5001f), -10);
219 TEST("rnd_halfintup(9.4999) == 9 ", vnl_math::rnd_halfintup(9.4999), 9);
220 TEST("rnd_halfintup(9.4999f) == 9 ", vnl_math::rnd_halfintup(9.4999f), 9);
221 TEST("rnd_halfintup(9.50) == 10", vnl_math::rnd_halfintup(9.50), 10);
222 TEST("rnd_halfintup(9.50f) == 10", vnl_math::rnd_halfintup(9.50f), 10);
223 TEST("rnd_halfintup(9.5001) == 10", vnl_math::rnd_halfintup(9.5001), 10);
224 TEST("rnd_halfintup(9.5001f) == 10", vnl_math::rnd_halfintup(9.5001f), 10);
225
226 TEST("floor(8.0) == 8", vnl_math::floor(8.0), 8);
227 TEST("floor(8.0f) == 8", vnl_math::floor(8.0f), 8);
228 TEST("floor(8.9999) == 8", vnl_math::floor(8.9999), 8);
229 TEST("floor(8.9999f) == 8", vnl_math::floor(8.9999f), 8);
230 TEST("floor(8.0001) == 8", vnl_math::floor(8.0001), 8);
231 TEST("floor(8.0001f) == 8", vnl_math::floor(8.0001f), 8);
232 TEST("floor(-8.0) == -8", vnl_math::floor(-8.0), -8);
233 TEST("floor(-8.0f) == -8", vnl_math::floor(-8.0f), -8);
234 TEST("floor(-8.9999) == -9", vnl_math::floor(-8.9999), -9);
235 TEST("floor(-8.9999f) == -9", vnl_math::floor(-8.9999f), -9);
236 TEST("floor(-8.0001) == -9", vnl_math::floor(-8.0001), -9);
237 TEST("floor(-8.0001f) == -9", vnl_math::floor(-8.0001f), -9);
238
239 TEST("floor(9.0) == 9 ", vnl_math::floor(9.0), 9);
240 TEST("floor(9.0f) == 9 ", vnl_math::floor(9.0f), 9);
241 TEST("floor(9.9999) == 9 ", vnl_math::floor(9.9999), 9);
242 TEST("floor(9.9999f) == 9 ", vnl_math::floor(9.9999f), 9);
243 TEST("floor(9.0001) == 9 ", vnl_math::floor(9.0001), 9);
244 TEST("floor(9.0001f) == 9 ", vnl_math::floor(9.0001f), 9);
245 TEST("floor(-9.0) == -9 ", vnl_math::floor(-9.0), -9);
246 TEST("floor(-9.0f) == -9 ", vnl_math::floor(-9.0f), -9);
247 TEST("floor(-9.9999) == -10", vnl_math::floor(-9.9999), -10);
248 TEST("floor(-9.9999f) == -10", vnl_math::floor(-9.9999f), -10);
249 TEST("floor(-9.0001) == -10", vnl_math::floor(-9.0001), -10);
250 TEST("floor(-9.0001f) == -10", vnl_math::floor(-9.0001f), -10);
251
252 TEST("ceil(8.0) == 8", vnl_math::ceil(8.0), 8);
253 TEST("ceil(8.0f) == 8", vnl_math::ceil(8.0f), 8);
254 TEST("ceil(8.9999) == 9", vnl_math::ceil(8.9999), 9);
255 TEST("ceil(8.9999f) == 9", vnl_math::ceil(8.9999f), 9);
256 TEST("ceil(8.0001) == 9", vnl_math::ceil(8.0001), 9);
257 TEST("ceil(8.0001f) == 9", vnl_math::ceil(8.0001f), 9);
258 TEST("ceil(-8.0) == -8", vnl_math::ceil(-8.0), -8);
259 TEST("ceil(-8.0f) == -8", vnl_math::ceil(-8.0f), -8);
260 TEST("ceil(-8.9999) == -8", vnl_math::ceil(-8.9999), -8);
261 TEST("ceil(-8.9999f) == -8", vnl_math::ceil(-8.9999f), -8);
262 TEST("ceil(-8.0001) == -8", vnl_math::ceil(-8.0001), -8);
263 TEST("ceil(-8.0001f) == -8", vnl_math::ceil(-8.0001f), -8);
264
265 TEST("ceil(9.0) == 9", vnl_math::ceil(9.0), 9);
266 TEST("ceil(9.0f) == 9", vnl_math::ceil(9.0f), 9);
267 TEST("ceil(9.9999) == 10", vnl_math::ceil(9.9999), 10);
268 TEST("ceil(9.9999f) == 10", vnl_math::ceil(9.9999f), 10);
269 TEST("ceil(9.0001) == 10", vnl_math::ceil(9.0001), 10);
270 TEST("ceil(9.0001f) == 10", vnl_math::ceil(9.0001f), 10);
271 TEST("ceil(-9.0) == -9", vnl_math::ceil(-9.0), -9);
272 TEST("ceil(-9.0f) == -9", vnl_math::ceil(-9.0f), -9);
273 TEST("ceil(-9.9999) == -9", vnl_math::ceil(-9.9999), -9);
274 TEST("ceil(-9.9999f) == -9", vnl_math::ceil(-9.9999f), -9);
275 TEST("ceil(-9.0001) == -9", vnl_math::ceil(-9.0001), -9);
276 TEST("ceil(-9.0001f) == -9", vnl_math::ceil(-9.0001f), -9);
277
278 TEST(" isfinite(f) ", vnl_math::isfinite(f), true);
279 TEST(" isfinite(d) ", vnl_math::isfinite(d), true);
280 TEST(" isfinite(i) ", vnl_math::isfinite(i), true);
281 TEST(" isfinite(z) ", vnl_math::isfinite(z), true);
282
283
284 // There is an assumption in this code that std::numeric_limits<float/double>::has_infinity==true
285
286 TEST("std::numeric_limits<float>::has_infinity==true assumption", std::numeric_limits<float>::has_infinity, true);
287 TEST("std::numeric_limits<double>::has_infinity==true assumption", std::numeric_limits<double>::has_infinity, true);
288 #ifdef INCLUDE_LONG_DOUBLE_TESTS
289 TEST("std::numeric_limits<ldouble>::has_infinity==true assumption",
290 std::numeric_limits<long double>::has_infinity,
291 true);
292 #endif
293 if (!std::numeric_limits<float>::has_infinity && !std::numeric_limits<double>::has_infinity)
294 {
295 std::cout << "Your platform doesn't appear to have an infinity. VXL is in places relatively\n"
296 << "dependent on the existence of an infinity. There are two solutions.\n"
297 << "A. If your platform really doesn't have an infinity, VXL's configuration code\n"
298 << " can be modified to correctly detect and use the infinity.\n"
299 << "B. Fix VXL so that it can cope with the lack of an infinity.\n"
300 << std::endl;
301 }
302 TEST("std::numeric_limits<float>::has_quiet_NaN==true assumption", std::numeric_limits<float>::has_quiet_NaN, true);
303 TEST("std::numeric_limits<double>::has_quiet_NaN==true assumption", std::numeric_limits<double>::has_quiet_NaN, true);
304 #ifdef INCLUDE_LONG_DOUBLE_TESTS
305 TEST("std::numeric_limits<ldouble>::has_quiet_NaN==true assumption",
306 std::numeric_limits<long double>::has_quiet_NaN,
307 true);
308 #endif
309 if (!std::numeric_limits<float>::has_quiet_NaN && !std::numeric_limits<double>::has_quiet_NaN)
310 {
311 std::cout << "Your platform doesn't appear to have a quiet NaN. VXL is in places relatively\n"
312 << "dependent on the existence of a quiet NaN. There are two solutions.\n"
313 << "A. If your platform really doesn't have a quiet NaN, VXL's configuration code\n"
314 << " can be modified to correctly detect and use the NaN.\n"
315 << "B. Fix VXL so that it can cope with the lack of a quiet NaN.\n"
316 << std::endl;
317 }
318 // Create Inf and -Inf:
319 const float pinf_f = std::numeric_limits<float>::infinity();
320 const float ninf_f = -std::numeric_limits<float>::infinity();
321 const double pinf_d = std::numeric_limits<double>::infinity();
322 const double ninf_d = -std::numeric_limits<double>::infinity();
323 // Create NaN
324 const float qnan_f = std::numeric_limits<float>::quiet_NaN();
325 const double qnan_d = std::numeric_limits<double>::quiet_NaN();
326 #ifdef INCLUDE_LONG_DOUBLE_TESTS
327 const long double pinf_q = std::numeric_limits<long double>::infinity();
328 const long double ninf_q = -std::numeric_limits<long double>::infinity();
329 const long double qnan_q = std::numeric_limits<long double>::quiet_NaN();
330 #endif
331
332 std::cout << "pinf_f = " << pinf_f << " = " << print_hex(pinf_f) << " sizeof " << sizeof(pinf_f) << '\n'
333 << "ninf_f = " << ninf_f << " = " << print_hex(ninf_f) << " sizeof " << sizeof(ninf_f) << '\n'
334 << "pinf_d = " << pinf_d << " = " << print_hex(pinf_d) << " sizeof " << sizeof(pinf_d) << '\n'
335 << "ninf_d = " << ninf_d << " = " << print_hex(ninf_d) << " sizeof " << sizeof(ninf_d) << '\n'
336 << "qnan_f = " << qnan_f << " = " << print_hex(qnan_f) << " sizeof " << sizeof(qnan_f) << '\n'
337 << "qnan_d = " << qnan_d << " = " << print_hex(qnan_d) << " sizeof " << sizeof(qnan_d) << '\n'
338 #ifdef INCLUDE_LONG_DOUBLE_TESTS
339 << "pinf_q = " << pinf_q << " = " << print_hex(pinf_q) << " sizeof " << sizeof(pinf_q) << '\n'
340 << "ninf_q = " << ninf_q << " = " << print_hex(ninf_q) << " sizeof " << sizeof(ninf_q) << '\n'
341 << "qnan_q = " << qnan_q << " = " << print_hex(qnan_q) << " sizeof " << sizeof(qnan_q) << '\n'
342 #endif
343 << std::endl;
344
345 TEST("!isfinite(pinf_f)", vnl_math::isfinite(pinf_f), false);
346 TEST("!isfinite(ninf_f)", vnl_math::isfinite(ninf_f), false);
347 TEST(" isinf(pinf_f) ", vnl_math::isinf(pinf_f), true);
348 TEST(" isinf(ninf_f) ", vnl_math::isinf(ninf_f), true);
349 TEST("!isnan(pinf_f) ", vnl_math::isnan(pinf_f), false);
350 TEST("!isnan(ninf_f) ", vnl_math::isnan(ninf_f), false);
351 TEST("!isfinite(qnan_f)", vnl_math::isfinite(qnan_f), false);
352 TEST("!isinf(qnan_f) ", vnl_math::isinf(qnan_f), false);
353 TEST(" isnan(qnan_f) ", vnl_math::isnan(qnan_f), true);
354
355 TEST("!isfinite(pinf_d)", vnl_math::isfinite(pinf_d), false);
356 TEST("!isfinite(ninf_d)", vnl_math::isfinite(ninf_d), false);
357 TEST(" isinf(pinf_d) ", vnl_math::isinf(pinf_d), true);
358 TEST(" isinf(ninf_d) ", vnl_math::isinf(ninf_d), true);
359 TEST("!isnan(pinf_d) ", vnl_math::isnan(pinf_d), false);
360 TEST("!isnan(ninf_d) ", vnl_math::isnan(ninf_d), false);
361 TEST("!isfinite(qnan_d)", vnl_math::isfinite(qnan_d), false);
362 TEST("!isinf(qnan_d) ", vnl_math::isinf(qnan_d), false);
363 TEST(" isnan(qnan_d) ", vnl_math::isnan(qnan_d), true);
364
365 #ifdef INCLUDE_LONG_DOUBLE_TESTS
366 # ifndef __ICC // "long double" has no standard internal representation on different platforms/compilers
367 TEST("!isfinite(pinf_q)", vnl_math::isfinite(pinf_q), false);
368 TEST("!isfinite(ninf_q)", vnl_math::isfinite(ninf_q), false);
369 TEST(" isinf(pinf_q) ", vnl_math::isinf(pinf_q), true);
370 TEST(" isinf(ninf_q) ", vnl_math::isinf(ninf_q), true);
371 TEST("!isnan(pinf_q) ", vnl_math::isnan(pinf_q), false);
372 TEST("!isnan(ninf_q) ", vnl_math::isnan(ninf_q), false);
373 TEST("!isfinite(qnan_q)", vnl_math::isfinite(qnan_q), false);
374 TEST("!isinf(qnan_q) ", vnl_math::isinf(qnan_q), false);
375 # endif // __ICC
376 #endif
377
378 TEST("!isfinite(huge_val(double))", vnl_math::isfinite(vnl_huge_val(double())), false);
379 TEST("!isfinite(huge_val(float))", vnl_math::isfinite(vnl_huge_val(float())), false);
380
381 // Test for math_sgn
382 TEST("vnl_math::sgn(+7) ", vnl_math::sgn(+7), 1);
383 TEST("vnl_math::sgn(-7) ", vnl_math::sgn(-7), -1);
384 TEST("vnl_math::sgn( 0) ", vnl_math::sgn(0), 0);
385
386 TEST("vnl_math::sgn(+7.0) ", vnl_math::sgn(+7.0), 1);
387 TEST("vnl_math::sgn(-7.0) ", vnl_math::sgn(-7.0), -1);
388 TEST("vnl_math::sgn(-0.0) ", vnl_math::sgn(-0.0), 0);
389 TEST("vnl_math::sgn(+0.0) ", vnl_math::sgn(-0.0), 0);
390
391 TEST("vnl_math::sgn(+7.0F) ", vnl_math::sgn(+7.0F), 1);
392 TEST("vnl_math::sgn(-7.0F) ", vnl_math::sgn(-7.0F), -1);
393 TEST("vnl_math::sgn(-0.0F) ", vnl_math::sgn(-0.0F), 0);
394 TEST("vnl_math::sgn(+0.0F) ", vnl_math::sgn(-0.0F), 0);
395
396 std::cout << std::endl;
397
398 // test vnl_math::angle_0_to_2pi() for "extreme values":
399 TEST("vnl_math::angle_0_to_2pi(2pi)", vnl_math::angle_0_to_2pi(vnl_math::twopi), 0.0);
400 double eps = vnl_math_test_2_epsilon; // which is smaller than the precision of vnl_math::pi
401 double conv_eps = vnl_math::angle_0_to_2pi(-eps);
402 std::cout << "conv_eps = " << conv_eps << " = 2pi - " << vnl_math::twopi - conv_eps << std::endl;
403 TEST("vnl_math::angle_0_to_2pi(-eps)", conv_eps < vnl_math::twopi && conv_eps > 6.283, true);
404 eps = vnl_math_test_20_epsilon; // which is larger than the precision of vnl_math::pi
405 conv_eps = vnl_math::angle_0_to_2pi(-eps);
406 std::cout << "conv_eps = " << conv_eps << " = 2pi - " << vnl_math::twopi - conv_eps << std::endl;
407 TEST(
408 "vnl_math::angle_0_to_2pi(-10eps)", conv_eps < vnl_math::twopi - vnl_math_test_2_epsilon && conv_eps > 6.283, true);
409 double ang = vnl_math::twopi - eps;
410 double conv_ang = vnl_math::angle_0_to_2pi(ang);
411 std::cout << "conv_ang = " << conv_ang << " = 2pi - " << vnl_math::twopi - conv_ang << std::endl;
412 TEST("vnl_math::angle_0_to_2pi(2pi-10eps)", conv_ang, ang);
413 // test vnl_math::angle_minuspi_to_pi() for "extreme values":
414 TEST("vnl_math::angle_minuspi_to_pi(2pi)", vnl_math::angle_minuspi_to_pi(vnl_math::twopi), 0.0);
415 TEST("vnl_math::angle_minuspi_to_pi(pi)", vnl_math::angle_minuspi_to_pi(vnl_math::pi), vnl_math::pi);
416 TEST("vnl_math::angle_minuspi_to_pi(-pi)", vnl_math::angle_minuspi_to_pi(-vnl_math::pi), -vnl_math::pi);
417 eps = vnl_math_test_2_epsilon; // which is smaller than the precision of vnl_math::pi
418 conv_eps = vnl_math::angle_minuspi_to_pi(-eps);
419 std::cout << "conv_eps = " << conv_eps << std::endl;
420 TEST("vnl_math::angle_minuspi_to_pi(-eps)", conv_eps, -eps);
421 eps = vnl_math_test_20_epsilon; // which is larger than the precision of vnl_math::pi
422 conv_eps = vnl_math::angle_minuspi_to_pi(-eps);
423 std::cout << "conv_eps = " << conv_eps << std::endl;
424 TEST("vnl_math::angle_minuspi_to_pi(-10eps)", conv_eps, -eps);
425
426 ///////////////
427 // TRUNCATED //
428 ///////////////
429
430 std::cout << "Truncated Remainder:" << std::endl;
431 // Truncated Remainder (% for integers, fmod for floating point)
432 // This behavior is most familiar to c++ programmers, but is unusual
433 // in mathematical terms.
434
435 {
436
437 std::cout << "+ +" << std::endl;
438
439 unsigned short x_short_u = 7;
440 unsigned short y_short_u = 2;
441 signed short x_short_s = 7;
442 signed short y_short_s = 2;
443 unsigned int x_int_u = 7;
444 unsigned int y_int_u = 2;
445 signed int x_int_s = 7;
446 signed int y_int_s = 2;
447 unsigned long x_long_u = 7;
448 unsigned long y_long_u = 2;
449 signed long x_long_s = 7;
450 signed long y_long_s = 2;
451 float x_float = 7;
452 float y_float = 2;
453 double x_double = 7;
454 double y_double = 2;
455 #ifdef INCLUDE_LONG_DOUBLE_TESTS
456 long double x_long_double = 7;
457 long double y_long_double = 2;
458 #endif
459
460 TEST("vnl_math::remainder_truncated(x_short_u ,y_short_u )",
461 vnl_math::remainder_truncated(x_short_u, y_short_u),
462 +1);
463 TEST("vnl_math::remainder_truncated(x_short_s ,y_short_s )",
464 vnl_math::remainder_truncated(x_short_s, y_short_s),
465 +1);
466 TEST("vnl_math::remainder_truncated(x_int_u ,y_int_u )",
467 vnl_math::remainder_truncated(x_int_u, y_int_u),
468 +1);
469 TEST("vnl_math::remainder_truncated(x_int_s ,y_int_s )",
470 vnl_math::remainder_truncated(x_int_s, y_int_s),
471 +1);
472 TEST("vnl_math::remainder_truncated(x_long_u ,y_long_u )",
473 vnl_math::remainder_truncated(x_long_u, y_long_u),
474 +1);
475 TEST("vnl_math::remainder_truncated(x_long_s ,y_long_s )",
476 vnl_math::remainder_truncated(x_long_s, y_long_s),
477 +1);
478 TEST("vnl_math::remainder_truncated(x_float ,y_float )",
479 vnl_math::remainder_truncated(x_float, y_float),
480 +1);
481 TEST("vnl_math::remainder_truncated(x_double ,y_double )",
482 vnl_math::remainder_truncated(x_double, y_double),
483 +1);
484 #ifdef INCLUDE_LONG_DOUBLE_TESTS
485 TEST("vnl_math::remainder_truncated(x_long_double,y_long_double)",
486 vnl_math::remainder_truncated(x_long_double, y_long_double),
487 +1);
488 #endif
489
490 std::cout << "+ -" << std::endl;
491
492 y_short_s *= -1;
493 y_int_s *= -1;
494 y_long_s *= -1;
495 y_float *= -1;
496 y_double *= -1;
497 #ifdef INCLUDE_LONG_DOUBLE_TESTS
498 y_long_double *= -1;
499 #endif
500
501 TEST("vnl_math::remainder_truncated(x_short_s ,y_short_s )",
502 vnl_math::remainder_truncated(x_short_s, y_short_s),
503 +1);
504 TEST("vnl_math::remainder_truncated(x_int_s ,y_int_s )",
505 vnl_math::remainder_truncated(x_int_s, y_int_s),
506 +1);
507 TEST("vnl_math::remainder_truncated(x_long_s ,y_long_s )",
508 vnl_math::remainder_truncated(x_long_s, y_long_s),
509 +1);
510 TEST("vnl_math::remainder_truncated(x_float ,y_float )",
511 vnl_math::remainder_truncated(x_float, y_float),
512 +1);
513 TEST("vnl_math::remainder_truncated(x_double ,y_double )",
514 vnl_math::remainder_truncated(x_double, y_double),
515 +1);
516 #ifdef INCLUDE_LONG_DOUBLE_TESTS
517 TEST("vnl_math::remainder_truncated(x_long_double,y_long_double)",
518 vnl_math::remainder_truncated(x_long_double, y_long_double),
519 +1);
520 #endif
521
522 std::cout << "- -" << std::endl;
523
524 x_short_s *= -1;
525 x_int_s *= -1;
526 x_long_s *= -1;
527 x_float *= -1;
528 x_double *= -1;
529 #ifdef INCLUDE_LONG_DOUBLE_TESTS
530 x_long_double *= -1;
531 #endif
532
533 TEST("vnl_math::remainder_truncated(x_short_s ,y_short_s )",
534 vnl_math::remainder_truncated(x_short_s, y_short_s),
535 -1);
536 TEST("vnl_math::remainder_truncated(x_int_s ,y_int_s )",
537 vnl_math::remainder_truncated(x_int_s, y_int_s),
538 -1);
539 TEST("vnl_math::remainder_truncated(x_long_s ,y_long_s )",
540 vnl_math::remainder_truncated(x_long_s, y_long_s),
541 -1);
542 TEST("vnl_math::remainder_truncated(x_float ,y_float )",
543 vnl_math::remainder_truncated(x_float, y_float),
544 -1);
545 TEST("vnl_math::remainder_truncated(x_double ,y_double )",
546 vnl_math::remainder_truncated(x_double, y_double),
547 -1);
548 #ifdef INCLUDE_LONG_DOUBLE_TESTS
549 TEST("vnl_math::remainder_truncated(x_long_double,y_long_double)",
550 vnl_math::remainder_truncated(x_long_double, y_long_double),
551 -1);
552 #endif
553
554 std::cout << "- +" << std::endl;
555
556 y_short_s *= -1;
557 y_int_s *= -1;
558 y_long_s *= -1;
559 y_float *= -1;
560 y_double *= -1;
561 #ifdef INCLUDE_LONG_DOUBLE_TESTS
562 y_long_double *= -1;
563 #endif
564
565 TEST("vnl_math::remainder_truncated(x_short_s ,y_short_s )",
566 vnl_math::remainder_truncated(x_short_s, y_short_s),
567 -1);
568 TEST("vnl_math::remainder_truncated(x_int_s ,y_int_s )",
569 vnl_math::remainder_truncated(x_int_s, y_int_s),
570 -1);
571 TEST("vnl_math::remainder_truncated(x_long_s ,y_long_s )",
572 vnl_math::remainder_truncated(x_long_s, y_long_s),
573 -1);
574 TEST("vnl_math::remainder_truncated(x_float ,y_float )",
575 vnl_math::remainder_truncated(x_float, y_float),
576 -1);
577 TEST("vnl_math::remainder_truncated(x_double ,y_double )",
578 vnl_math::remainder_truncated(x_double, y_double),
579 -1);
580 #ifdef INCLUDE_LONG_DOUBLE_TESTS
581 TEST("vnl_math::remainder_truncated(x_long_double,y_long_double)",
582 vnl_math::remainder_truncated(x_long_double, y_long_double),
583 -1);
584 #endif
585 }
586
587 /////////////
588 // FLOORED //
589 /////////////
590
591 std::cout << "Floored Remainder:" << std::endl;
592 // This definition is identical to truncated_remainder for two positive arguments.
593 // When one or both of the arguments are negative, this attempts to mimick Python's modulus behavior.
594
595 {
596
597 std::cout << "+ +" << std::endl;
598
599 unsigned short x_short_u = 7;
600 unsigned short y_short_u = 2;
601 signed short x_short_s = 7;
602 signed short y_short_s = 2;
603 unsigned int x_int_u = 7;
604 unsigned int y_int_u = 2;
605 signed int x_int_s = 7;
606 signed int y_int_s = 2;
607 unsigned long x_long_u = 7;
608 unsigned long y_long_u = 2;
609 signed long x_long_s = 7;
610 signed long y_long_s = 2;
611 float x_float = 7;
612 float y_float = 2;
613 double x_double = 7;
614 double y_double = 2;
615 #ifdef INCLUDE_LONG_DOUBLE_TESTS
616 long double x_long_double = 7;
617 long double y_long_double = 2;
618 #endif
619
620 TEST("vnl_math::remainder_floored(x_short_u ,y_short_u )",
621 vnl_math::remainder_floored(x_short_u, y_short_u),
622 +1);
623 TEST("vnl_math::remainder_floored(x_short_s ,y_short_s )",
624 vnl_math::remainder_floored(x_short_s, y_short_s),
625 +1);
626 TEST("vnl_math::remainder_floored(x_int_u ,y_int_u )", vnl_math::remainder_floored(x_int_u, y_int_u), +1);
627 TEST("vnl_math::remainder_floored(x_int_s ,y_int_s )", vnl_math::remainder_floored(x_int_s, y_int_s), +1);
628 TEST(
629 "vnl_math::remainder_floored(x_long_u ,y_long_u )", vnl_math::remainder_floored(x_long_u, y_long_u), +1);
630 TEST(
631 "vnl_math::remainder_floored(x_long_s ,y_long_s )", vnl_math::remainder_floored(x_long_s, y_long_s), +1);
632 TEST("vnl_math::remainder_floored(x_float ,y_float )", vnl_math::remainder_floored(x_float, y_float), +1);
633 TEST(
634 "vnl_math::remainder_floored(x_double ,y_double )", vnl_math::remainder_floored(x_double, y_double), +1);
635 #ifdef INCLUDE_LONG_DOUBLE_TESTS
636 TEST("vnl_math::remainder_floored(x_long_double,y_long_double)",
637 vnl_math::remainder_floored(x_long_double, y_long_double),
638 +1);
639 #endif
640
641 std::cout << "+ -" << std::endl;
642
643 y_short_s *= -1;
644 y_int_s *= -1;
645 y_long_s *= -1;
646 y_float *= -1;
647 y_double *= -1;
648 #ifdef INCLUDE_LONG_DOUBLE_TESTS
649 y_long_double *= -1;
650 #endif
651
652 TEST("vnl_math::remainder_floored(x_short_s ,y_short_s )",
653 vnl_math::remainder_floored(x_short_s, y_short_s),
654 -1);
655 TEST("vnl_math::remainder_floored(x_int_s ,y_int_s )", vnl_math::remainder_floored(x_int_s, y_int_s), -1);
656 TEST(
657 "vnl_math::remainder_floored(x_long_s ,y_long_s )", vnl_math::remainder_floored(x_long_s, y_long_s), -1);
658 TEST("vnl_math::remainder_floored(x_float ,y_float )", vnl_math::remainder_floored(x_float, y_float), -1);
659 TEST(
660 "vnl_math::remainder_floored(x_double ,y_double )", vnl_math::remainder_floored(x_double, y_double), -1);
661 #ifdef INCLUDE_LONG_DOUBLE_TESTS
662 TEST("vnl_math::remainder_floored(x_long_double,y_long_double)",
663 vnl_math::remainder_floored(x_long_double, y_long_double),
664 -1);
665 #endif
666
667 std::cout << "- -" << std::endl;
668
669 x_short_s *= -1;
670 x_int_s *= -1;
671 x_long_s *= -1;
672 x_float *= -1;
673 x_double *= -1;
674 #ifdef INCLUDE_LONG_DOUBLE_TESTS
675 x_long_double *= -1;
676 #endif
677
678 TEST("vnl_math::remainder_floored(x_short_s ,y_short_s )",
679 vnl_math::remainder_floored(x_short_s, y_short_s),
680 -1);
681 TEST("vnl_math::remainder_floored(x_int_s ,y_int_s )", vnl_math::remainder_floored(x_int_s, y_int_s), -1);
682 TEST(
683 "vnl_math::remainder_floored(x_long_s ,y_long_s )", vnl_math::remainder_floored(x_long_s, y_long_s), -1);
684 TEST("vnl_math::remainder_floored(x_float ,y_float )", vnl_math::remainder_floored(x_float, y_float), -1);
685 TEST(
686 "vnl_math::remainder_floored(x_double ,y_double )", vnl_math::remainder_floored(x_double, y_double), -1);
687 #ifdef INCLUDE_LONG_DOUBLE_TESTS
688 TEST("vnl_math::remainder_floored(x_long_double,y_long_double)",
689 vnl_math::remainder_floored(x_long_double, y_long_double),
690 -1);
691 #endif
692
693 std::cout << "- +" << std::endl;
694
695 y_short_s *= -1;
696 y_int_s *= -1;
697 y_long_s *= -1;
698 y_float *= -1;
699 y_double *= -1;
700 #ifdef INCLUDE_LONG_DOUBLE_TESTS
701 y_long_double *= -1;
702 #endif
703
704 TEST("vnl_math::remainder_floored(x_short_s ,y_short_s )",
705 vnl_math::remainder_floored(x_short_s, y_short_s),
706 +1);
707 TEST("vnl_math::remainder_floored(x_int_s ,y_int_s )", vnl_math::remainder_floored(x_int_s, y_int_s), +1);
708 TEST(
709 "vnl_math::remainder_floored(x_long_s ,y_long_s )", vnl_math::remainder_floored(x_long_s, y_long_s), +1);
710 TEST("vnl_math::remainder_floored(x_float ,y_float )", vnl_math::remainder_floored(x_float, y_float), +1);
711 TEST(
712 "vnl_math::remainder_floored(x_double ,y_double )", vnl_math::remainder_floored(x_double, y_double), +1);
713 #ifdef INCLUDE_LONG_DOUBLE_TESTS
714 TEST("vnl_math::remainder_floored(x_long_double,y_long_double)",
715 vnl_math::remainder_floored(x_long_double, y_long_double),
716 +1);
717 #endif
718 }
719
720 #define RETURN_TYPE_TEST(funcname, argtypename, returntypename) \
721 { \
722 const bool test_return_type = \
723 std::is_same<decltype(vnl_math ::funcname(static_cast<argtypename>(123.4))), returntypename>(); \
724 TEST("vnl_math::" #funcname "<" #argtypename "> returns " #returntypename " type", test_return_type, true); \
725 } \
726 void()
727
728 RETURN_TYPE_TEST(isinf, int, bool);
729 RETURN_TYPE_TEST(isinf, float, bool);
730 RETURN_TYPE_TEST(isinf, double, bool);
731
732 RETURN_TYPE_TEST(isnan, int, bool);
733 RETURN_TYPE_TEST(isnan, float, bool);
734 RETURN_TYPE_TEST(isnan, double, bool);
735
736 RETURN_TYPE_TEST(isfinite, int, bool);
737 RETURN_TYPE_TEST(isfinite, float, bool);
738 RETURN_TYPE_TEST(isfinite, double, bool);
739
740 RETURN_TYPE_TEST(isnormal, int, bool);
741 RETURN_TYPE_TEST(isnormal, float, bool);
742 RETURN_TYPE_TEST(isnormal, double, bool);
743 }
744
745 TESTMAIN(test_math);
746