1 #include <iostream>
2
3 #include "testlib/testlib_test.h"
4 #include "vnl/vnl_real_npolynomial.h"
5
6 void
test_real_npolynomial()7 test_real_npolynomial()
8 {
9 vnl_vector<double> coef_0(3), coef_1(3), coef_2(4);
10 for (unsigned int i = 0; i < 3; ++i)
11 coef_0(i) = i + 1.0; // f0 = X + 2Y + 3Z
12 for (unsigned int i = 0; i < 3; ++i)
13 coef_1(i) = i + 1.0; // f1 = X^2 + 2XY^3 + 3
14 for (unsigned int i = 0; i < 4; ++i)
15 coef_2(i) = 2 * i + 1.0; // f2 = X^3 + 3X^2Y + 5XY^2 + 7Y^3
16
17 vnl_matrix<unsigned int> expo_0(3, 3, 0U), expo_1(3, 2, 0U), expo_2(4, 2, 0U);
18 for (unsigned int i = 0; i < 3; ++i)
19 expo_0(i, i) = 1;
20 expo_1(0, 0) = 2;
21 expo_1(1, 0) = 1;
22 expo_1(1, 1) = 3;
23 for (unsigned int i = 0; i < 4; ++i)
24 expo_2(i, 1) = expo_2(3 - i, 0) = i;
25
26 vnl_real_npolynomial f0(coef_0, expo_0), f1(coef_1, expo_1), f2(coef_2, expo_2);
27
28 std::cout << "f0 = " << f0 << "f1 = " << f1 << "f2 = " << f2;
29 TEST("f0 has total degree 1", f0.degree(), 1);
30 TEST("f0 has maximal degree 1", f0.maxdegree(), 1);
31 TEST("f0 has degree 1 in X", f0.degrees()[0], 1);
32 TEST("f0 has degree 1 in Y", f0.degrees()[1], 1);
33 TEST("f0 has degree 1 in Z", f0.degrees()[2], 1);
34
35 TEST("f1 has total degree 4", f1.degree(), 4);
36 TEST("f1 has maximal degree 3", f1.maxdegree(), 3);
37 TEST("f1 has degree 2 in X", f1.degrees()[0], 2);
38 TEST("f1 has degree 3 in Y", f1.degrees()[1], 3);
39
40 TEST("f2 has total degree 3", f2.degree(), 3);
41 TEST("f2 has maximal degree 3", f2.maxdegree(), 3);
42 TEST("f2 has degree 3 in X", f2.degrees()[0], 3);
43 TEST("f2 has degree 3 in Y", f2.degrees()[1], 3);
44
45 vnl_vector<double> vec3(2);
46 vec3(0) = vec3(1) = 2.5;
47
48 vnl_real_npolynomial f3 = f1 + f2;
49 std::cout << "f1+f2 = " << f3;
50 TEST("f1+f2=f3", f1.eval(vec3) + f2.eval(vec3), f3.eval(vec3));
51
52 TEST("f3 has total degree 4", f3.degree(), 4);
53 TEST("f3 has maximal degree 3", f3.maxdegree(), 3);
54 TEST("f3 has degree 3 in X", f3.degrees()[0], 3);
55 TEST("f3 has degree 3 in Y", f3.degrees()[1], 3);
56
57 vnl_real_npolynomial f4 = f1 - f2;
58 std::cout << "f1-f2 = " << f4;
59 TEST("f1-f2=f4", f1.eval(vec3) - f2.eval(vec3), f4.eval(vec3));
60
61 TEST("f4 has total degree 4", f4.degree(), 4);
62 TEST("f4 has maximal degree 3", f4.maxdegree(), 3);
63 TEST("f4 has degree 3 in X", f4.degrees()[0], 3);
64 TEST("f4 has degree 3 in Y", f4.degrees()[1], 3);
65
66 vnl_real_npolynomial f5 = f1 * f2;
67 std::cout << "f1*f2 = " << f5;
68 TEST("f1*f2=f5", f1.eval(vec3) * f2.eval(vec3), f5.eval(vec3));
69
70 TEST("f5 has total degree 7", f5.degree(), 7);
71 TEST("f5 has maximal degree 6", f5.maxdegree(), 6);
72 TEST("f5 has degree 5 in X", f5.degrees()[0], 5);
73 TEST("f5 has degree 6 in Y", f5.degrees()[1], 6);
74
75 TEST("f1*f2 has correct total degree", f5.degree(), f1.degree() + f2.degree());
76 TEST("f1*f2 has correct maximal degree", f5.maxdegree(), f1.maxdegree() + f2.maxdegree());
77 TEST("f1*f2 has correct degree in X", f5.degrees()[0], f1.degrees()[0] + f2.degrees()[0]);
78 TEST("f1*f2 has correct degree in Y", f5.degrees()[1], f1.degrees()[1] + f2.degrees()[1]);
79 }
80
81 TESTMAIN(test_real_npolynomial);
82