1 /*
2 Copyright (C) 2011 Fredrik Johansson
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 <stdio.h>
13 #include <stdlib.h>
14 #include <gmp.h>
15 #include "flint.h"
16 #include "arith.h"
17 #include "fmpz.h"
18 #include "fmpz_poly.h"
19 #include "ulong_extras.h"
20
cyclotomic_naive(fmpz_poly_t poly,ulong n)21 void cyclotomic_naive(fmpz_poly_t poly, ulong n)
22 {
23 fmpz_poly_t t;
24 slong d;
25
26 fmpz_poly_init(t);
27
28 fmpz_poly_set_ui(poly, UWORD(1));
29 for (d = 1; d <= n; d++)
30 {
31 if (n % d == 0)
32 {
33 if (n_moebius_mu(n / d) == 1)
34 {
35 fmpz_poly_zero(t);
36 fmpz_poly_set_coeff_si(t, d, 1);
37 fmpz_poly_set_coeff_si(t, 0, -1);
38 fmpz_poly_mul(poly, poly, t);
39 }
40 }
41 }
42
43 for (d = 1; d <= n; d++)
44 {
45 if (n % d == 0)
46 {
47 if (n_moebius_mu(n / d) == -1)
48 {
49 fmpz_poly_zero(t);
50 fmpz_poly_set_coeff_si(t, d, 1);
51 fmpz_poly_set_coeff_si(t, 0, -1);
52 fmpz_poly_div(poly, poly, t);
53 }
54 }
55 }
56
57 fmpz_poly_clear(t);
58 }
59
main()60 int main()
61 {
62 fmpz_poly_t A, B;
63 slong n;
64
65 FLINT_TEST_INIT(state);
66
67 flint_printf("cyclotomic_polynomial....");
68 fflush(stdout);
69
70 for (n = 0; n <= 1000; n++)
71 {
72 fmpz_poly_init(A);
73 fmpz_poly_init(B);
74
75 arith_cyclotomic_polynomial(A, n);
76 cyclotomic_naive(B, n);
77
78 if (!fmpz_poly_equal(A, B))
79 {
80 flint_printf("FAIL: wrong value of Phi_%wd(x)\n", n);
81 flint_printf("Computed:\n");
82 fmpz_poly_print_pretty(A, "x");
83 flint_printf("\n\nExpected:\n");
84 fmpz_poly_print_pretty(B, "x");
85 flint_printf("\n\n");
86 abort();
87 }
88
89 fmpz_poly_clear(A);
90 fmpz_poly_clear(B);
91 }
92
93 /* We verify the first value that does not fit on 32 bits.
94 This exercises the slow path at least on a 32 bit system.
95 Testing the 64 bit value is a bit too much to do by default
96 as it requires ~2 GB of memory and takes a few minutes. */
97 {
98 fmpz_t h, ref;
99
100 const ulong nn = UWORD(10163195);
101 /* const ulong nn = UWORD(169828113); 64-bit case */
102
103 fmpz_init(h);
104 fmpz_init(ref);
105 fmpz_set_str(ref, "1376877780831", 10);
106 /* fmpz_set_str(ref, "31484567640915734941", 10); 64-bit case */
107
108 fmpz_poly_init(A);
109 arith_cyclotomic_polynomial(A, UWORD(10163195));
110 fmpz_poly_height(h, A);
111
112 if (!fmpz_equal(h, ref))
113 {
114 flint_printf("Bad computation of Phi_%wd(x)\n", nn);
115 flint_printf("Computed height:\n");
116 fmpz_print(h);
117 flint_printf("\nExpected height:\n");
118 fmpz_print(ref);
119 flint_printf("\n\n");
120 abort();
121 }
122
123 fmpz_poly_clear(A);
124 fmpz_clear(h);
125 fmpz_clear(ref);
126 }
127
128 FLINT_TEST_CLEANUP(state);
129 flint_printf("PASS\n");
130 return 0;
131 }
132