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