1 /*
2 Copyright (C) 2016 Vincent Delecroix
3
4 This file is part of FLINT.
5
6 FLINT is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include <gmp.h>
13 #include "flint.h"
14 #include "fmpz.h"
15 #include "fmpz_poly.h"
16
17 int
main(void)18 main(void)
19 {
20 int i, j, k, l, result;
21 fmpz_t il, jl, kl, tot;
22 slong n;
23 fmpz_poly_t a, b, c, d, e, f;
24
25 FLINT_TEST_INIT(state);
26
27 flint_printf("power_sums....");
28
29 /* Check that it is valid in degree 3 with integer roots, ie */
30 /* for polynomials of the form (x-i)(x-j)(x-k) */
31 for (i = -5; i < 5; i++)
32 for (j = i; j < 5; j++)
33 for (k = j; k < 5; k++)
34 {
35 fmpz_poly_init(a);
36 fmpz_poly_init(b);
37 fmpz_poly_init(c);
38
39 fmpz_poly_set_coeff_si(a, 0, -i * j * k);
40 fmpz_poly_set_coeff_si(a, 1, i * j + i * k + j * k);
41 fmpz_poly_set_coeff_si(a, 2, -i - j - k);
42 fmpz_poly_set_coeff_si(a, 3, 1);
43
44 fmpz_poly_power_sums(b, a, 20);
45 fmpz_poly_power_sums_naive(c, a, 20);
46 result = fmpz_poly_equal(b, c);
47 if (!result)
48 {
49 flint_printf("FAIL:\n");
50 flint_printf("%d %d %d\n\n", i, j, k);
51 fmpz_poly_print(b), flint_printf("\n\n");
52 fmpz_poly_print(c), flint_printf("\n\n");
53 abort();
54 }
55
56 fmpz_init(il);
57 fmpz_init(jl);
58 fmpz_init(kl);
59 fmpz_set_ui(il, 1);
60 fmpz_set_ui(jl, 1);
61 fmpz_set_ui(kl, 1);
62 fmpz_init(tot);
63 for (l = 0; l < FLINT_MIN(20, fmpz_poly_length(b)); l++)
64 {
65 fmpz_set(tot, il);
66 fmpz_add(tot, tot, jl);
67 fmpz_add(tot, tot, kl);
68 result = fmpz_equal(fmpz_poly_get_coeff_ptr(b, l), tot);
69 if (!result)
70 {
71 flint_printf("FAIL:\n\n");
72 flint_printf("%d %d %d %d\n", i, j, k, l);
73 abort();
74 }
75 fmpz_mul_si(il, il, i);
76 fmpz_mul_si(jl, jl, j);
77 fmpz_mul_si(kl, kl, k);
78 }
79
80 fmpz_poly_power_sums(b, a, 4);
81 fmpz_poly_power_sums_to_poly(c, b);
82 result = fmpz_poly_equal(a, c);
83 if (!result)
84 {
85 flint_printf("FAIL: newton series to poly\n\n");
86 fmpz_poly_print(a), flint_printf("\n\n");
87 fmpz_poly_print(b), flint_printf("\n\n");
88 fmpz_poly_print(c), flint_printf("\n\n");
89 }
90
91 fmpz_clear(il);
92 fmpz_clear(jl);
93 fmpz_clear(kl);
94 fmpz_clear(tot);
95 fmpz_poly_clear(a);
96 fmpz_poly_clear(b);
97 fmpz_poly_clear(c);
98 }
99
100 /* Check that the various implementations coincide and that */
101 /* power_sums_to_poly gives back the original polynomial */
102 for (i = 0; i < 20 * flint_test_multiplier(); i++)
103 {
104 fmpz_poly_init(a);
105 fmpz_poly_init(b);
106 fmpz_poly_init(c);
107 fmpz_poly_init(d);
108 fmpz_poly_init(e);
109 fmpz_poly_init(f);
110 fmpz_poly_randtest_not_zero(a, state, 1 + n_randint(state, 50), 100);
111 fmpz_poly_set_coeff_ui(a, fmpz_poly_degree(a), 1);
112
113 for (n = -1; n < 3 * fmpz_poly_degree(a);
114 n += 1 + n_randint(state, fmpz_poly_degree(a)))
115 {
116 fmpz_poly_power_sums(b, a, n);
117 fmpz_poly_power_sums_naive(c, a, n);
118
119 result = fmpz_poly_equal(b, c);
120 if (!result)
121 {
122 flint_printf
123 ("FAIL: equality power_sums, power_sums_naive\n");
124 flint_printf("%ld", n);
125 fmpz_poly_print(a), flint_printf("\n\n");
126 abort();
127 }
128
129 if (n >= 1)
130 {
131 /* use the formula PowerSums(p) = rev(poly') / rev(poly) */
132 fmpz_poly_reverse(e, a, fmpz_poly_length(a));
133 fmpz_poly_derivative(f, a);
134 fmpz_poly_reverse(f, f, fmpz_poly_length(f));
135 fmpz_poly_div_series(d, f, e, n);
136 result = fmpz_poly_equal(b, d);
137 if (!result)
138 {
139 flint_printf("FAIL: equality with schoenhage formula\n");
140 flint_printf("%ld", n);
141 fmpz_poly_print(a), flint_printf("\n\n");
142 abort();
143 }
144 }
145 }
146
147 fmpz_poly_power_sums(b, a, fmpz_poly_length(a));
148 fmpz_poly_power_sums_to_poly(c, b);
149 result = fmpz_poly_equal(a, c);
150 if (!result)
151 {
152 flint_printf("FAIL: power_sums_to_poly\n");
153 fmpz_poly_print(a), flint_printf("\n\n");
154 fmpz_poly_print(b), flint_printf("\n\n");
155 fmpz_poly_print(c), flint_printf("\n\n");
156 }
157
158 fmpz_poly_clear(a);
159 fmpz_poly_clear(b);
160 fmpz_poly_clear(c);
161 fmpz_poly_clear(d);
162 fmpz_poly_clear(e);
163 fmpz_poly_clear(f);
164 }
165
166 /* Check that the product of polynomials correspond to the sum of Power sums series */
167 for (i = 0; i < 20 * flint_test_multiplier(); i++)
168 {
169 fmpz_poly_init(a);
170 fmpz_poly_init(b);
171 fmpz_poly_init(c);
172 fmpz_poly_init(d);
173 fmpz_poly_init(e);
174 fmpz_poly_init(f);
175
176 fmpz_poly_randtest_not_zero(a, state, 1 + n_randint(state, 5), 20);
177 fmpz_poly_set_coeff_ui(a, fmpz_poly_degree(a), 1);
178 fmpz_poly_randtest_not_zero(b, state, 1 + n_randint(state, 5), 20);
179 fmpz_poly_set_coeff_ui(b, fmpz_poly_degree(b), 1);
180
181 fmpz_poly_power_sums(c, a, 20);
182 fmpz_poly_power_sums(d, b, 20);
183 fmpz_poly_add(f, c, d);
184
185 fmpz_poly_mul(e, a, b);
186 fmpz_poly_power_sums(e, e, 20);
187
188 result = fmpz_poly_equal(e, f);
189 if (!result)
190 {
191 flint_printf("FAIL: PowerSums(p1 p2) = PowerSums(p1) + PowerSums(p2)\n");
192 fmpz_poly_print(a), flint_printf("\n");
193 fmpz_poly_print(b), flint_printf("\n");
194 fmpz_poly_print(c), flint_printf("\n");
195 fmpz_poly_print(d), flint_printf("\n");
196 fmpz_poly_print(e), flint_printf("\n");
197 fmpz_poly_print(f), flint_printf("\n");
198 abort();
199 }
200
201 fmpz_poly_clear(a);
202 fmpz_poly_clear(b);
203 fmpz_poly_clear(c);
204 fmpz_poly_clear(d);
205 fmpz_poly_clear(e);
206 fmpz_poly_clear(f);
207 }
208
209 FLINT_TEST_CLEANUP(state);
210
211 flint_printf("PASS\n");
212 return 0;
213 }
214