1 #include <check.h>
2 #include <mps/mps.h>
3
START_TEST(test_chebyshev_poly_80)4 START_TEST (test_chebyshev_poly_80)
5 {
6 mps_context * ctx = mps_context_new ();
7 mps_chebyshev_poly *cp = mps_chebyshev_poly_new (ctx, 80, MPS_STRUCTURE_REAL_INTEGER);
8 mpc_t *mroots = NULL;
9 rdpe_t *radii = NULL;
10 mpq_t one, zero;
11 int i;
12
13 mpq_init (one);
14 mpq_init (zero);
15
16 mpq_set_ui (one, 1U, 1U);
17 mpq_set_ui (zero, 0U, 1U);
18
19 mps_chebyshev_poly_set_coefficient_q (ctx, cp, 80, one, zero);
20 for (i = 0; i < 80; i++)
21 {
22 mps_chebyshev_poly_set_coefficient_q (ctx, cp, i, zero, zero);
23 }
24
25 mps_context_set_input_poly (ctx, MPS_POLYNOMIAL (cp));
26 mps_context_select_algorithm (ctx, MPS_ALGORITHM_SECULAR_GA);
27 mps_thread_pool_set_concurrency_limit (ctx, NULL, 1);
28 mps_mpsolve (ctx);
29
30 /* Check that the roots are correct */
31 mps_context_get_roots_m (ctx, &mroots, &radii);
32 for (i = 0; i < 80; i++)
33 {
34 int j, found_root = -1;
35 rdpe_t expected_root;
36 rdpe_t diff;
37 double epsilon = DBL_MAX;
38
39 rdpe_set_d (expected_root, cos ((2.0 * i + 1) / 160 * PI));;
40
41 printf ("[Chebyshev tests] Expected root %d: (%1.20lf, 0)\n",
42 i, rdpe_get_d (expected_root));
43
44 for (j = 0; j < 80; j++)
45 {
46 cdpe_t ctmp;
47 double residue;
48
49 mpc_get_cdpe (ctmp, mroots[j]);
50 rdpe_sub (diff, expected_root, cdpe_Re (ctmp));
51
52 residue = sqrt (pow (fabs (rdpe_get_d (diff)), 2) + pow (fabs (rdpe_get_d (cdpe_Im (ctmp))), 2));
53 if (residue < epsilon)
54 {
55 epsilon = residue;
56 found_root = j;
57 }
58 }
59
60 printf ("[Chebyshev tests] Residue for approximation %3d: %e\n", i, epsilon);
61 printf ("[Chebyshev tests] Inclusion radii: %e\n", rdpe_get_d (radii[found_root]));
62 fail_unless (epsilon < 4.0 * DBL_EPSILON + rdpe_get_d (radii[found_root]));
63 }
64
65 mps_polynomial_free (ctx, MPS_POLYNOMIAL (cp));
66 mps_context_free (ctx);
67
68 mpq_clear (one);
69 mpq_clear (zero);
70 }
71 END_TEST
72
START_TEST(test_chebyshev_poly_20)73 START_TEST (test_chebyshev_poly_20)
74 {
75 mps_context * ctx = mps_context_new ();
76 mps_chebyshev_poly *cp = mps_chebyshev_poly_new (ctx, 20, MPS_STRUCTURE_REAL_INTEGER);
77 mpc_t *mroots = NULL;
78 rdpe_t *radii = NULL;
79 mpq_t one, zero;
80 int i;
81
82 mpq_init (one);
83 mpq_init (zero);
84
85 mpq_set_ui (one, 1U, 1U);
86 mpq_set_ui (zero, 0U, 1U);
87
88 mps_chebyshev_poly_set_coefficient_q (ctx, cp, 20, one, zero);
89 for (i = 0; i < 20; i++)
90 {
91 mps_chebyshev_poly_set_coefficient_q (ctx, cp, i, zero, zero);
92 }
93
94 mps_context_set_input_poly (ctx, MPS_POLYNOMIAL (cp));
95 mps_context_select_algorithm (ctx, MPS_ALGORITHM_SECULAR_GA);
96 mps_thread_pool_set_concurrency_limit (ctx, NULL, 1);
97 mps_mpsolve (ctx);
98
99 /* Check that the roots are correct */
100 mps_context_get_roots_m (ctx, &mroots, &radii);
101 for (i = 0; i < 20; i++)
102 {
103 int j, found_root = -1;
104 rdpe_t expected_root;
105 rdpe_t diff;
106 double epsilon = DBL_MAX;
107
108 rdpe_set_d (expected_root, cos ((2.0 * i + 1) / 40 * PI));;
109
110 printf ("[Chebyshev tests] Expected root %d: (%1.20lf, 0)\n",
111 i, rdpe_get_d (expected_root));
112
113 for (j = 0; j < 20; j++)
114 {
115 cdpe_t ctmp;
116 double residue;
117
118 mpc_get_cdpe (ctmp, mroots[j]);
119 rdpe_sub (diff, expected_root, cdpe_Re (ctmp));
120
121 residue = sqrt (pow (fabs (rdpe_get_d (diff)), 2) + pow (fabs (rdpe_get_d (cdpe_Im (ctmp))), 2));
122 if (residue < epsilon)
123 {
124 epsilon = residue;
125 found_root = j;
126 }
127 }
128
129 printf ("[Chebyshev tests] Residue for approximation %3d: %e\n", i, epsilon);
130 printf ("[Chebyshev tests] Inclusion radii: %e\n", rdpe_get_d (radii[found_root]));
131 fail_unless (epsilon < 4.0 * DBL_EPSILON + rdpe_get_d (radii[found_root]));
132 }
133
134 mps_polynomial_free (ctx, MPS_POLYNOMIAL (cp));
135 mps_context_free (ctx);
136
137 mpq_clear (one);
138 mpq_clear (zero);
139 }
140 END_TEST
141
142 Suite*
chebyshev_suite(void)143 chebyshev_suite (void)
144 {
145 Suite *s = suite_create ("chebyshev");
146
147 TCase *tcase_t = tcase_create ("Solution of Chebyshev polynomials");
148
149 tcase_add_test (tcase_t, test_chebyshev_poly_20);
150 tcase_add_test (tcase_t, test_chebyshev_poly_80);
151
152 suite_add_tcase (s, tcase_t);
153 return s;
154 }
155
156 int
main(void)157 main (void)
158 {
159 Suite *cs = chebyshev_suite ();
160 SRunner *sr = srunner_create (cs);
161 int number_failed;
162
163 if (getenv ("MPS_SIMPLE_TESTS_ONLY") != NULL)
164 return EXIT_SUCCESS;
165
166 srunner_run_all (sr, CK_NORMAL);
167
168 /* Get number of failed test and report */
169 number_failed = srunner_ntests_failed (sr);
170 srunner_free (sr);
171
172 return(number_failed != 0);
173 }
174