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 <math.h>
13 #include "arith.h"
14 
15 #define MAX_32BIT 58
16 
17 static const int lookup_table[MAX_32BIT][28] =
18 {
19  {-1, 1}, {1, 1}, {1, 2}, {0, 1}, {-1, 2, 4}, {-1, 2},
20  {-1, -4, 4, 8}, {-1, 0, 2}, {1, -6, 0, 8}, {-1, -2, 4},
21  {1, 6, -12, -32, 16, 32}, {-3, 0, 4}, {-1, 6, 24, -32, -80, 32, 64},
22  {1, -4, -4, 8}, {1, 8, -16, -8, 16}, {1, 0, -8, 0, 8},
23  {1, -8, -40, 80, 240, -192, -448, 128, 256}, {-1, -6, 0, 8},
24  {1, 10, -40, -160, 240, 672, -448, -1024, 256, 512}, {5, 0, -20, 0, 16},
25  {1, -16, 32, 48, -96, -32, 64}, {-1, 6, 12, -32, -16, 32},
26  {-1, -12, 60, 280, -560, -1792, 1792, 4608, -2304, -5120, 1024, 2048},
27  {1, 0, -16, 0, 16}, {-1, 10, 100, -40, -800, 32, 2240, 0, -2560, 0,
28   1024}, {-1, -6, 24, 32, -80, -32, 64},
29  {1, 18, 0, -240, 0, 864, 0, -1152, 0, 512}, {-7, 0, 56, 0, -112, 0, 64},
30  {-1, 14, 112, -448, -2016, 4032, 13440, -15360, -42240, 28160, 67584,
31   -24576, -53248, 8192, 16384}, {1, -8, -16, 8, 16},
32  {-1, -16, 112, 672, -2016, -8064, 13440, 42240, -42240, -112640, 67584,
33   159744, -53248, -114688, 16384, 32768},
34  {1, 0, -32, 0, 160, 0, -256, 0, 128},
35  {1, -24, 48, 344, -688, -1088, 2176, 1280, -2560, -512, 1024},
36  {1, 8, -40, -80, 240, 192, -448, -128, 256},
37  {1, 16, -160, -368, 1760, 2272, -7232, -5504, 13824, 5632, -12288,
38   -2048, 4096}, {-3, 0, 36, 0, -96, 0, 64},
39  {-1, 18, 180, -960, -5280, 14784, 59136, -101376, -329472, 366080,
40   1025024, -745472, -1863680, 860160, 1966080, -524288, -1114112, 131072,
41   262144}, {-1, 10, 40, -160, -240, 672, 448, -1024, -256, 512},
42  {1, 24, -48, -632, 1264, 3296, -6592, -6784, 13568, 6144, -12288, -2048,
43   4096}, {1, 0, -48, 0, 304, 0, -512, 0, 256},
44  {1, -20, -220, 1320, 7920, -25344, -109824, 219648, 768768, -1025024,
45   -3075072, 2795520, 7454720, -4587520, -11141120, 4456448, 10027008,
46   -2359296, -4980736, 524288, 1048576}, {1, 16, 32, -48, -96, 32, 64},
47  {1, 22, -220, -1760, 7920, 41184, -109824, -439296, 768768, 2562560,
48   -3075072, -8945664, 7454720, 19496960, -11141120, -26738688, 10027008,
49   22413312, -4980736, -10485760, 1048576, 2097152},
50  {-11, 0, 220, 0, -1232, 0, 2816, 0, -2816, 0, 1024},
51  {1, -24, -144, 248, 1680, -864, -7168, 1152, 13824, -512, -12288, 0,
52   4096}, {1, -12, -60, 280, 560, -1792, -1792, 4608, 2304, -5120, -1024,
53   2048}, {-1, -24, 264, 2288, -11440, -64064, 192192, 823680, -1647360,
54   -5857280, 8200192, 25346048, -25346048, -70189056, 50135040, 127008768,
55   -63504384, -149422080, 49807360, 110100480, -22020096, -46137344,
56   4194304, 8388608}, {1, 0, -64, 0, 320, 0, -512, 0, 256},
57  {-1, 28, 196, -2968, -3136, 66304, 18816, -658816, -53760, 3587584,
58   78848, -11741184, -57344, 24084480, 16384, -31195136, 0, 24772608, 0,
59   -11010048, 0, 2097152}, {-1, -10, 100, 40, -800, -32, 2240, 0, -2560,
60   0, 1024}, {1, 32, -64, -1504, 3008, 16832, -33664, -76288, 152576,
61   173568, -347136, -210944, 421888, 131072, -262144, -32768, 65536},
62  {13, 0, -364, 0, 2912, 0, -9984, 0, 16640, 0, -13312, 0, 4096},
63  {-1, 26, 364, -2912, -21840, 96096, 512512, -1464320, -6223360,
64   12446720, 44808192, -65175552, -206389248, 222265344, 635043840,
65   -508035072, -1333592064, 784465920, 1917583360, -807403520,
66   -1857028096, 530579456, 1157627904, -201326592, -419430400, 33554432,
67   67108864}, {-1, 18, 0, -240, 0, 864, 0, -1152, 0, 512},
68  {1, 24, -432, -1208, 15216, 28064, -185024, -263424, 1149184, 1250304,
69   -4177920, -3356672, 9375744, 5324800, -13123584, -4947968, 11141120,
70   2490368, -5242880, -524288, 1048576},
71  {1, 0, -96, 0, 1376, 0, -6656, 0, 13568, 0, -12288, 0, 4096},
72  {1, -40, 80, 2120, -4240, -31648, 63296, 194432, -388864, -613376,
73   1226752, 1087488, -2174976, -1097728, 2195456, 589824, -1179648,
74   -131072, 262144}, {-1, -14, 112, 448, -2016, -4032, 13440, 15360,
75   -42240, -28160, 67584, 24576, -53248, -8192, 16384}
76 };
77 
78 /* The coefficients in 2^d * \prod_{i=1}^d (x - cos(a_i)) are
79    easily bounded using the binomial theorem. */
80 static slong
magnitude_bound(slong d)81 magnitude_bound(slong d)
82 {
83     slong res;
84     fmpz_t t;
85     fmpz_init(t);
86     fmpz_bin_uiui(t, d, d / 2);
87     res = fmpz_bits(t);
88     fmpz_clear(t);
89     return FLINT_ABS(res) + d;
90 }
91 
92 static void
fmpz_mul_or_div_2exp(fmpz_t x,fmpz_t y,slong s)93 fmpz_mul_or_div_2exp(fmpz_t x, fmpz_t y, slong s)
94 {
95     if (s >= 0)
96         fmpz_mul_2exp(x, y, s);
97     else
98         fmpz_fdiv_q_2exp(x, y, -s);
99 }
100 
101 
102 /* Balanced product of linear factors (x+alpha_i) using
103    fixed-point arithmetic with prec bits */
104 static void
balanced_product(fmpz * c,fmpz * alpha,slong len,slong prec)105 balanced_product(fmpz * c, fmpz * alpha, slong len, slong prec)
106 {
107     if (len == 1)
108     {
109         fmpz_one(c + 1);
110         fmpz_mul_2exp(c + 1, c + 1, prec);
111         fmpz_set(c, alpha);
112     }
113     else if (len == 2)
114     {
115         fmpz_mul(c, alpha, alpha + 1);
116         fmpz_fdiv_q_2exp(c, c, prec);
117         fmpz_add(c + 1, alpha, alpha + 1);
118         fmpz_one(c + 2);
119         fmpz_mul_2exp(c + 2, c + 2, prec);
120     }
121     else
122     {
123         fmpz *L, *R;
124         slong i, m;
125 
126         m = len / 2;
127         L = _fmpz_vec_init(len + 2);
128         R = L + m + 1;
129 
130         balanced_product(L, alpha, m, prec);
131         balanced_product(R, alpha + m, len - m, prec);
132         _fmpz_poly_mul(c, R, len - m + 1, L, m + 1);
133 
134         for (i = 0; i < len + 1; i++)
135             fmpz_fdiv_q_2exp(c + i, c + i, prec);
136 
137         _fmpz_vec_clear(L, len + 2);
138     }
139 }
140 
141 void
_arith_cos_minpoly(fmpz * coeffs,slong d,ulong n)142 _arith_cos_minpoly(fmpz * coeffs, slong d, ulong n)
143 {
144     slong i, j;
145     fmpz * alpha;
146     fmpz_t half;
147     mpfr_t t, u;
148     flint_bitcnt_t prec;
149     slong exp;
150 
151     if (n <= MAX_32BIT)
152     {
153         for (i = 0; i <= d; i++)
154             fmpz_set_si(coeffs + i, lookup_table[n - 1][i]);
155         return;
156     }
157 
158     /* Direct formula for odd primes > 3 */
159     if (n_is_prime(n))
160     {
161         slong s = (n - 1) / 2;
162 
163         switch (s % 4)
164         {
165             case 0:
166                 fmpz_set_si(coeffs, WORD(1));
167                 fmpz_set_si(coeffs + 1, -s);
168                 break;
169             case 1:
170                 fmpz_set_si(coeffs, WORD(1));
171                 fmpz_set_si(coeffs + 1, s + 1);
172                 break;
173             case 2:
174                 fmpz_set_si(coeffs, WORD(-1));
175                 fmpz_set_si(coeffs + 1, s);
176                 break;
177             case 3:
178                 fmpz_set_si(coeffs, WORD(-1));
179                 fmpz_set_si(coeffs + 1, -s - 1);
180                 break;
181         }
182 
183         for (i = 2; i <= s; i++)
184         {
185             slong b = (s - i) % 2;
186             fmpz_mul2_uiui(coeffs + i, coeffs + i - 2, s+i-b, s+2-b-i);
187             fmpz_divexact2_uiui(coeffs + i, coeffs + i, i, i-1);
188             fmpz_neg(coeffs + i, coeffs + i);
189         }
190 
191         return;
192     }
193 
194     prec = magnitude_bound(d) + 5 + FLINT_BIT_COUNT(d);
195 
196     alpha = _fmpz_vec_init(d);
197     fmpz_init(half);
198     mpfr_init2(t, prec);
199     mpfr_init2(u, prec);
200 
201     fmpz_one(half);
202     fmpz_mul_2exp(half, half, prec - 1);
203     mpfr_const_pi(t, prec);
204     mpfr_div_ui(t, t, n, MPFR_RNDN);
205 
206     for (i = j = 0; j < d; i++)
207     {
208         if (n_gcd(n, i) == 1)
209         {
210             mpfr_mul_ui(u, t, 2 * i, MPFR_RNDN);
211             mpfr_cos(u, u, MPFR_RNDN);
212             mpfr_neg(u, u, MPFR_RNDN);
213             exp = mpfr_get_z_2exp(_fmpz_promote(alpha + j), u);
214             _fmpz_demote_val(alpha + j);
215             fmpz_mul_or_div_2exp(alpha + j, alpha + j, exp + prec);
216             j++;
217         }
218     }
219 
220     balanced_product(coeffs, alpha, d, prec);
221 
222     /* Scale and round */
223     for (i = 0; i < d + 1; i++)
224     {
225         slong r = d;
226         if ((n & (n - 1)) == 0)
227             r--;
228         fmpz_mul_2exp(coeffs + i, coeffs + i, r);
229         fmpz_add(coeffs + i, coeffs + i, half);
230         fmpz_fdiv_q_2exp(coeffs + i, coeffs + i, prec);
231     }
232 
233     fmpz_clear(half);
234     mpfr_clear(t);
235     mpfr_clear(u);
236     _fmpz_vec_clear(alpha, d);
237 }
238 
239 void
arith_cos_minpoly(fmpz_poly_t poly,ulong n)240 arith_cos_minpoly(fmpz_poly_t poly, ulong n)
241 {
242     if (n == 0)
243     {
244         fmpz_poly_set_ui(poly, UWORD(1));
245     }
246     else
247     {
248         slong d = (n <= 2) ? 1 : n_euler_phi(n) / 2;
249 
250         fmpz_poly_fit_length(poly, d + 1);
251         _arith_cos_minpoly(poly->coeffs, d, n);
252         _fmpz_poly_set_length(poly, d + 1);
253     }
254 }
255