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 "fmpq.h"
16 #include "fmpz_poly.h"
17 #include "fmpq_poly.h"
18
19 int
main(void)20 main(void)
21 {
22 int i, j, k, l, den, result;
23 fmpz_t il, jl, kl, dl;
24 fmpq_t tot;
25 fmpq_t tmp;
26 fmpq_poly_t a, b, c, d;
27 fmpz_poly_t az, bz;
28
29 FLINT_TEST_INIT(state);
30
31 flint_printf("power_sums....");
32
33 /* Check that it is valid in degree 3 with rational roots, ie */
34 /* for polynomials of the form (dx-i)(dx-j)(dx-k) */
35 for (den = 1; den < 6; den++)
36 for (i = -4; i < 4; i++)
37 for (j = i; j < 4; j++)
38 for (k = j; k < 4; k++)
39 {
40 fmpq_poly_init(a);
41 fmpq_poly_init(b);
42 fmpq_poly_init(c);
43
44 fmpq_poly_set_coeff_si(a, 0, -i * j * k);
45 fmpq_poly_set_coeff_si(a, 1,
46 den * (i * j + i * k + j * k));
47 fmpq_poly_set_coeff_si(a, 2, den * den * (-i - j - k));
48 fmpq_poly_set_coeff_si(a, 3, den * den * den);
49
50 fmpq_poly_power_sums(b, a, 20);
51
52 fmpz_init(il);
53 fmpz_init(jl);
54 fmpz_init(kl);
55 fmpz_init(dl);
56 fmpq_init(tot);
57 fmpq_init(tmp);
58 fmpz_one(il);
59 fmpz_one(jl);
60 fmpz_one(kl);
61 fmpz_one(dl);
62 for (l = 0; l < FLINT_MIN(20, fmpq_poly_length(b)); l++)
63 {
64 fmpq_zero(tot);
65 fmpq_add_fmpz(tot, tot, il);
66 fmpq_add_fmpz(tot, tot, jl);
67 fmpq_add_fmpz(tot, tot, kl);
68 fmpq_div_fmpz(tot, tot, dl);
69
70 fmpq_poly_get_coeff_fmpq(tmp, b, l);
71 result = fmpq_equal(tmp, tot);
72 if (!result)
73 {
74 flint_printf("FAIL: power sums rational root\n");
75 flint_printf("%d %d %d %d %d\n", i, j, k, den, l);
76 abort();
77 }
78 fmpz_mul_si(il, il, i);
79 fmpz_mul_si(jl, jl, j);
80 fmpz_mul_si(kl, kl, k);
81 fmpz_mul_si(dl, dl, den);
82 }
83
84 fmpq_poly_power_sums(b, a, 4);
85 fmpq_poly_power_sums_to_poly(c, b);
86 fmpq_poly_make_monic(a, a);
87 result = fmpq_poly_equal(a, c);
88 if (!result)
89 {
90 flint_printf("FAIL: power sums to poly\n");
91 fmpq_poly_print(a), flint_printf("\n\n");
92 fmpq_poly_print(b), flint_printf("\n\n");
93 abort();
94 }
95
96 fmpz_clear(il);
97 fmpz_clear(jl);
98 fmpz_clear(kl);
99 fmpq_clear(tot);
100 fmpq_clear(tmp);
101 fmpq_poly_clear(a);
102 fmpq_poly_clear(b);
103 fmpq_poly_clear(c);
104 }
105
106 /* Check that going back and forth between the power sums representation gives the identity */
107 for (i = 0; i < 50 * flint_test_multiplier(); i++)
108 {
109 fmpz_poly_init(az);
110 fmpz_poly_init(bz);
111 fmpq_poly_init(a);
112 fmpq_poly_init(b);
113 fmpq_poly_init(c);
114
115 fmpz_poly_randtest_not_zero(az, state, 1 + n_randint(state, 15), 30);
116 fmpq_poly_set_fmpz_poly(a, az);
117
118 fmpq_poly_power_sums(b, a, 20);
119 fmpq_poly_power_sums_to_fmpz_poly(bz, b);
120 fmpq_poly_power_sums_to_poly(c, b);
121
122 fmpq_poly_make_monic(a, a);
123 _fmpz_poly_primitive_part(az->coeffs, az->coeffs, az->length);
124
125 result = fmpq_poly_equal(a, c) && fmpz_poly_equal(bz, az);
126 if (!result)
127 {
128 flint_printf("FAIL: power sums - power sums to poly\n");
129 fmpz_poly_print(az), flint_printf("\n\n");
130 fmpz_poly_print(bz), flint_printf("\n\n");
131 fmpq_poly_print(a), flint_printf("\n\n");
132 fmpq_poly_print(b), flint_printf("\n\n");
133 fmpq_poly_print(c), flint_printf("\n\n");
134 abort();
135 }
136
137 fmpq_poly_clear(a);
138 fmpq_poly_clear(b);
139 fmpq_poly_clear(c);
140 fmpz_poly_clear(az);
141 fmpz_poly_clear(bz);
142 }
143
144
145 /* Check that the product of polynomials correspond to the sum of Power sums series */
146 /* (and aliasing of fmpq_poly_power_sums) */
147 for (i = 0; i < 20 * flint_test_multiplier(); i++)
148 {
149 fmpq_poly_init(a);
150 fmpq_poly_init(b);
151 fmpq_poly_init(c);
152 fmpq_poly_init(d);
153
154 fmpq_poly_randtest_not_zero(a, state, 1 + n_randint(state, 10), 30);
155 fmpq_poly_randtest_not_zero(b, state, 1 + n_randint(state, 10), 30);
156
157 fmpq_poly_mul(c, a, b);
158 fmpq_poly_power_sums(c, c, 20);
159
160 fmpq_poly_power_sums(a, a, 20);
161 fmpq_poly_power_sums(b, b, 20);
162 fmpq_poly_add(d, a, b);
163
164 result = fmpq_poly_equal(c, d);
165 if (!result)
166 {
167 flint_printf
168 ("FAIL: PowerSums(p1 p2) = PowerSums(p1) + PowerSums(p2)\n");
169 fmpq_poly_print(a), flint_printf("\n");
170 fmpq_poly_print(b), flint_printf("\n");
171 fmpq_poly_print(c), flint_printf("\n");
172 fmpq_poly_print(d), flint_printf("\n");
173 abort();
174 }
175
176 fmpq_poly_clear(a);
177 fmpq_poly_clear(b);
178 fmpq_poly_clear(c);
179 fmpq_poly_clear(d);
180 }
181
182 FLINT_TEST_CLEANUP(state);
183
184 flint_printf("PASS\n");
185 return 0;
186 }
187