mirror of
https://github.com/vale981/arb
synced 2025-03-05 17:31:38 -05:00
arb-based cosine minpoly calculation
This commit is contained in:
Fredrik Johansson
parent
39fd03828f
commit
9c6e1d71f3
6 changed files with 472 additions and 3 deletions
|
|
@ -84,3 +84,30 @@ The exponential function
|
|||
This implementation guarantees consistent rounding but will overflow
|
||||
for too large *x*.
|
||||
|
||||
Trigonometric functions
|
||||
--------------------------------------------------------------------------------
|
||||
|
||||
.. function:: void _elefun_cos_minpoly_roots(fmprb_struct * alpha, long d, ulong n, long prec)
|
||||
|
||||
Sets the vector *alpha* to the the *d* roots of `\Phi_n(x)`, computed
|
||||
using a working precision of *prec* bits.
|
||||
|
||||
.. function:: void _elefun_cos_minpoly(fmpz * coeffs, long d, ulong n)
|
||||
|
||||
.. function:: void elefun_cos_minpoly(fmpz_poly_t poly, ulong n)
|
||||
|
||||
Computes `\Phi_n(x)`, the minimal polynomial of `\cos(2\pi/n)`,
|
||||
which for `n > 2` has degree `d = \varphi(n) / 2`.
|
||||
For small `n`, the coefficients of the polynomial are looked up
|
||||
from a table. For large `n`, we compute numerical approximations of
|
||||
the roots using :func:`_elefun_cos_minpoly_roots`, then multiply
|
||||
the linear factors together using a balanced product tree, and convert
|
||||
the numerical coefficients to exact integers.
|
||||
|
||||
Since `\Phi_n(x) = 2^r \prod_{i=1}^d (x - \cos(\alpha_i))` for some
|
||||
`\alpha_i`, where `r = d - 1` if `n` is a power of two and `r = d`
|
||||
otherwise, we can use the binomial theorem to estimate the required
|
||||
working precision as `d + \log_2 {d \choose d / 2}`, plus a few
|
||||
guard bits. This estimate is not optimal, but it is acceptably tight
|
||||
for large `n`.
|
||||
|
||||
|
|
|
|||
8
elefun.h
8
elefun.h
|
|
@ -27,6 +27,7 @@
|
|||
#define ELEFUN_H
|
||||
|
||||
#include "fmprb.h"
|
||||
#include "fmprb_poly.h"
|
||||
#include "fmpz_extras.h"
|
||||
|
||||
static __inline__ void
|
||||
|
|
@ -57,5 +58,12 @@ int elefun_exp_precomp(fmprb_t z, const fmprb_t x, long prec, int minus_one);
|
|||
|
||||
void elefun_exp_via_mpfr(fmprb_t z, const fmprb_t x, long prec);
|
||||
|
||||
|
||||
void _elefun_cos_minpoly_roots(fmprb_struct * alpha, long d, ulong n, long prec);
|
||||
|
||||
void _elefun_cos_minpoly(fmpz * coeffs, long d, ulong n);
|
||||
|
||||
void elefun_cos_minpoly(fmpz_poly_t poly, ulong n);
|
||||
|
||||
#endif
|
||||
|
||||
|
|
|
|||
188
elefun/cos_minpoly.c
Normal file
188
elefun/cos_minpoly.c
Normal file
|
|
@ -0,0 +1,188 @@
|
|||
/*=============================================================================
|
||||
|
||||
This file is part of ARB.
|
||||
|
||||
ARB is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
ARB is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with ARB; if not, write to the Free Software
|
||||
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
|
||||
|
||||
=============================================================================*/
|
||||
/******************************************************************************
|
||||
|
||||
Copyright (C) 2013 Fredrik Johansson
|
||||
|
||||
******************************************************************************/
|
||||
|
||||
#include "elefun.h"
|
||||
#include "ulong_extras.h"
|
||||
|
||||
#define ELEFUN_COS_MINPOLY_MAX_LOOKUP 58
|
||||
#define MAX_LENGTH 28
|
||||
|
||||
static const int lookup_table[ELEFUN_COS_MINPOLY_MAX_LOOKUP][MAX_LENGTH] =
|
||||
{
|
||||
{-1, 1},
|
||||
{1, 1},
|
||||
{1, 2},
|
||||
{0, 1},
|
||||
{-1, 2, 4},
|
||||
{-1, 2},
|
||||
{-1, -4, 4, 8},
|
||||
{-1, 0, 2},
|
||||
{1, -6, 0, 8},
|
||||
{-1, -2, 4},
|
||||
{1, 6, -12, -32, 16, 32},
|
||||
{-3, 0, 4},
|
||||
{-1, 6, 24, -32, -80, 32, 64},
|
||||
{1, -4, -4, 8},
|
||||
{1, 8, -16, -8, 16},
|
||||
{1, 0, -8, 0, 8},
|
||||
{1, -8, -40, 80, 240, -192, -448, 128, 256},
|
||||
{-1, -6, 0, 8},
|
||||
{1, 10, -40, -160, 240, 672, -448, -1024, 256, 512},
|
||||
{5, 0, -20, 0, 16},
|
||||
{1, -16, 32, 48, -96, -32, 64},
|
||||
{-1, 6, 12, -32, -16, 32},
|
||||
{-1, -12, 60, 280, -560, -1792, 1792, 4608, -2304, -5120, 1024, 2048},
|
||||
{1, 0, -16, 0, 16},
|
||||
{-1, 10, 100, -40, -800, 32, 2240, 0, -2560, 0, 1024},
|
||||
{-1, -6, 24, 32, -80, -32, 64},
|
||||
{1, 18, 0, -240, 0, 864, 0, -1152, 0, 512},
|
||||
{-7, 0, 56, 0, -112, 0, 64},
|
||||
{-1, 14, 112, -448, -2016, 4032, 13440, -15360, -42240, 28160, 67584,
|
||||
-24576, -53248, 8192, 16384},
|
||||
{1, -8, -16, 8, 16},
|
||||
{-1, -16, 112, 672, -2016, -8064, 13440, 42240, -42240, -112640, 67584,
|
||||
159744, -53248, -114688, 16384, 32768},
|
||||
{1, 0, -32, 0, 160, 0, -256, 0, 128},
|
||||
{1, -24, 48, 344, -688, -1088, 2176, 1280, -2560, -512, 1024},
|
||||
{1, 8, -40, -80, 240, 192, -448, -128, 256},
|
||||
{1, 16, -160, -368, 1760, 2272, -7232, -5504, 13824, 5632, -12288,
|
||||
-2048, 4096},
|
||||
{-3, 0, 36, 0, -96, 0, 64},
|
||||
{-1, 18, 180, -960, -5280, 14784, 59136, -101376, -329472, 366080,
|
||||
1025024, -745472, -1863680, 860160, 1966080, -524288, -1114112, 131072,
|
||||
262144},
|
||||
{-1, 10, 40, -160, -240, 672, 448, -1024, -256, 512},
|
||||
{1, 24, -48, -632, 1264, 3296, -6592, -6784, 13568, 6144, -12288, -2048,
|
||||
4096},
|
||||
{1, 0, -48, 0, 304, 0, -512, 0, 256},
|
||||
{1, -20, -220, 1320, 7920, -25344, -109824, 219648, 768768, -1025024,
|
||||
-3075072, 2795520, 7454720, -4587520, -11141120, 4456448, 10027008,
|
||||
-2359296, -4980736, 524288, 1048576},
|
||||
{1, 16, 32, -48, -96, 32, 64},
|
||||
{1, 22, -220, -1760, 7920, 41184, -109824, -439296, 768768, 2562560,
|
||||
-3075072, -8945664, 7454720, 19496960, -11141120, -26738688, 10027008,
|
||||
22413312, -4980736, -10485760, 1048576, 2097152},
|
||||
{-11, 0, 220, 0, -1232, 0, 2816, 0, -2816, 0, 1024},
|
||||
{1, -24, -144, 248, 1680, -864, -7168, 1152, 13824, -512, -12288, 0, 4096},
|
||||
{1, -12, -60, 280, 560, -1792, -1792, 4608, 2304, -5120, -1024, 2048},
|
||||
{-1, -24, 264, 2288, -11440, -64064, 192192, 823680, -1647360,
|
||||
-5857280, 8200192, 25346048, -25346048, -70189056, 50135040, 127008768,
|
||||
-63504384, -149422080, 49807360, 110100480, -22020096, -46137344,
|
||||
4194304, 8388608},
|
||||
{1, 0, -64, 0, 320, 0, -512, 0, 256},
|
||||
{-1, 28, 196, -2968, -3136, 66304, 18816, -658816, -53760, 3587584, 78848,
|
||||
-11741184, -57344, 24084480, 16384, -31195136, 0, 24772608, 0, -11010048,
|
||||
0, 2097152},
|
||||
{-1, -10, 100, 40, -800, -32, 2240, 0, -2560, 0, 1024},
|
||||
{1, 32, -64, -1504, 3008, 16832, -33664, -76288, 152576, 173568, -347136,
|
||||
-210944, 421888, 131072, -262144, -32768, 65536},
|
||||
{13, 0, -364, 0, 2912, 0, -9984, 0, 16640, 0, -13312, 0, 4096},
|
||||
{-1, 26, 364, -2912, -21840, 96096, 512512, -1464320, -6223360, 12446720,
|
||||
44808192, -65175552, -206389248, 222265344, 635043840, -508035072,
|
||||
-1333592064, 784465920, 1917583360, -807403520, -1857028096, 530579456,
|
||||
1157627904, -201326592, -419430400, 33554432, 67108864},
|
||||
{-1, 18, 0, -240, 0, 864, 0, -1152, 0, 512},
|
||||
{1, 24, -432, -1208, 15216, 28064, -185024, -263424, 1149184, 1250304,
|
||||
-4177920, -3356672, 9375744, 5324800, -13123584, -4947968, 11141120,
|
||||
2490368, -5242880, -524288, 1048576},
|
||||
{1, 0, -96, 0, 1376, 0, -6656, 0, 13568, 0, -12288, 0, 4096},
|
||||
{1, -40, 80, 2120, -4240, -31648, 63296, 194432, -388864, -613376,
|
||||
1226752, 1087488, -2174976, -1097728, 2195456, 589824, -1179648,
|
||||
-131072, 262144}, {-1, -14, 112, 448, -2016, -4032, 13440, 15360,
|
||||
-42240, -28160, 67584, 24576, -53248, -8192, 16384}
|
||||
};
|
||||
|
||||
/* The coefficients in 2^d * \prod_{i=1}^d (x - cos(a_i)) are
|
||||
easily bounded using the binomial theorem. */
|
||||
static long
|
||||
magnitude_bound(long d)
|
||||
{
|
||||
long res;
|
||||
fmpz_t t;
|
||||
fmpz_init(t);
|
||||
fmpz_bin_uiui(t, d, d / 2);
|
||||
res = fmpz_bits(t);
|
||||
fmpz_clear(t);
|
||||
return FLINT_ABS(res) + d;
|
||||
}
|
||||
|
||||
void
|
||||
_elefun_cos_minpoly(fmpz * coeffs, long d, ulong n)
|
||||
{
|
||||
fmprb_struct * alpha, * fcoeffs;
|
||||
fmpz_t t;
|
||||
len_t i, prec;
|
||||
|
||||
if (n <= ELEFUN_COS_MINPOLY_MAX_LOOKUP)
|
||||
{
|
||||
for (i = 0; i <= d; i++)
|
||||
fmpz_set_si(coeffs + i, lookup_table[n - 1][i]);
|
||||
}
|
||||
else
|
||||
{
|
||||
fmpz_init(t);
|
||||
alpha = _fmprb_vec_init(d);
|
||||
fcoeffs = _fmprb_vec_init(d + 1);
|
||||
|
||||
prec = magnitude_bound(d) + 10 + FLINT_BIT_COUNT(d);
|
||||
|
||||
_elefun_cos_minpoly_roots(alpha, d, n, prec);
|
||||
_fmprb_poly_product_roots(fcoeffs, alpha, d, prec);
|
||||
|
||||
for (i = 0; i <= d; i++)
|
||||
{
|
||||
fmprb_mul_2exp_si(fcoeffs + i, fcoeffs + i,
|
||||
(n & (n - 1)) == 0 ? d - 1 : d);
|
||||
|
||||
if (!fmprb_get_unique_fmpz(coeffs + i, fcoeffs + i))
|
||||
{
|
||||
printf("insufficient precision computing cos minpolynomial!\n");
|
||||
abort();
|
||||
}
|
||||
}
|
||||
|
||||
_fmprb_vec_clear(alpha, d);
|
||||
_fmprb_vec_clear(fcoeffs, d + 1);
|
||||
fmpz_clear(t);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
elefun_cos_minpoly(fmpz_poly_t poly, ulong n)
|
||||
{
|
||||
if (n == 0)
|
||||
{
|
||||
fmpz_poly_set_si(poly, 1);
|
||||
}
|
||||
else
|
||||
{
|
||||
long d = (n <= 2) ? 1 : n_euler_phi(n) / 2;
|
||||
|
||||
fmpz_poly_fit_length(poly, d + 1);
|
||||
_elefun_cos_minpoly(poly->coeffs, d, n);
|
||||
_fmpz_poly_set_length(poly, d + 1);
|
||||
}
|
||||
}
|
||||
|
||||
77
elefun/cos_minpoly_roots.c
Normal file
77
elefun/cos_minpoly_roots.c
Normal file
|
|
@ -0,0 +1,77 @@
|
|||
/*=============================================================================
|
||||
|
||||
This file is part of ARB.
|
||||
|
||||
ARB is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
ARB is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with ARB; if not, write to the Free Software
|
||||
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
|
||||
|
||||
=============================================================================*/
|
||||
/******************************************************************************
|
||||
|
||||
Copyright (C) 2013 Fredrik Johansson
|
||||
|
||||
******************************************************************************/
|
||||
|
||||
#include "elefun.h"
|
||||
#include "ulong_extras.h"
|
||||
|
||||
/* computes all the d roots alpha_j of Phi_n(x) */
|
||||
void
|
||||
_elefun_cos_minpoly_roots(fmprb_struct * alpha, long d, ulong n, long prec)
|
||||
{
|
||||
fmprb_t t, u, v, s1, s2, c1, c2;
|
||||
long i, j;
|
||||
|
||||
fmprb_init(t);
|
||||
fmprb_init(u);
|
||||
fmprb_init(v);
|
||||
fmprb_init(s1);
|
||||
fmprb_init(s2);
|
||||
fmprb_init(c1);
|
||||
fmprb_init(c2);
|
||||
|
||||
fmprb_const_pi(s1, prec);
|
||||
fmprb_mul_ui(s1, s1, 2, prec);
|
||||
fmprb_div_ui(s1, s1, n, prec);
|
||||
fmprb_sin_cos(s1, c1, s1, prec);
|
||||
fmprb_one(c2);
|
||||
fmprb_zero(s2);
|
||||
|
||||
for (i = j = 0; j < d; i++)
|
||||
{
|
||||
if (n_gcd(n, i) == 1)
|
||||
{
|
||||
fmprb_set(alpha + j, c2);
|
||||
j++;
|
||||
}
|
||||
|
||||
fmprb_mul(t, c1, c2, prec);
|
||||
fmprb_mul(u, s1, s2, prec);
|
||||
fmprb_sub(v, t, u, prec);
|
||||
|
||||
fmprb_mul(t, c1, s2, prec);
|
||||
fmprb_mul(u, s1, c2, prec);
|
||||
fmprb_add(s2, t, u, prec);
|
||||
fmprb_set(c2, v);
|
||||
}
|
||||
|
||||
fmprb_clear(t);
|
||||
fmprb_clear(u);
|
||||
fmprb_clear(v);
|
||||
fmprb_clear(s1);
|
||||
fmprb_clear(s2);
|
||||
fmprb_clear(c1);
|
||||
fmprb_clear(c2);
|
||||
}
|
||||
|
||||
169
elefun/test/t-cos_minpoly.c
Normal file
169
elefun/test/t-cos_minpoly.c
Normal file
|
|
@ -0,0 +1,169 @@
|
|||
/*=============================================================================
|
||||
|
||||
This file is part of FLINT.
|
||||
|
||||
FLINT is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
FLINT is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with FLINT; if not, write to the Free Software
|
||||
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
|
||||
|
||||
=============================================================================*/
|
||||
/******************************************************************************
|
||||
|
||||
Copyright (C) 2011, 013 Fredrik Johansson
|
||||
|
||||
******************************************************************************/
|
||||
|
||||
#include "elefun.h"
|
||||
|
||||
/*
|
||||
Generated with Mathematica:
|
||||
Table[Mod[MinimalPolynomial[Cos[2 Pi/n]][1337], 31337], {n,1,500}]
|
||||
*/
|
||||
|
||||
static const short testdata[] = {
|
||||
1,
|
||||
1336, 1338, 2675, 1337, 8113, 2673, 6283, 2719, 29508, 2765, 6949,
|
||||
5437, 2788, 26742, 25554, 26194, 29376, 29506, 30945, 15614, 8957,
|
||||
16643, 9263, 21050, 30556, 10533, 1570, 11562, 3988, 16546, 26642, 4041,
|
||||
3581, 109, 9839, 27175, 11691, 1460, 28287, 18369, 16503, 3184, 13336,
|
||||
23083, 12495, 3246, 14160, 8081, 5301, 8652, 28989, 24149, 17733, 1568,
|
||||
4800, 28863, 29280, 13741, 30919, 29819, 28584, 8913, 550, 6207, 13930,
|
||||
23373, 12644, 15265, 27975, 30386, 1603, 15894, 22276, 3138, 11610,
|
||||
2208, 515, 30817, 23050, 4333, 25031, 13615, 5116, 18609, 25490, 14555,
|
||||
22663, 8425, 21751, 19293, 3, 10688, 26829, 14467, 1426, 12413, 5305,
|
||||
25377, 27164, 3711,
|
||||
9613, 22340, 7457, 3704, 1795, 22877, 31060, 17472, 11317, 22274,
|
||||
11036, 7796, 27242, 22174, 3663, 10507, 16599, 18192, 15208, 7257, 7022,
|
||||
10810, 27891, 18495, 7032, 11383, 20768, 27351, 31089, 27723, 10486,
|
||||
2075, 25298, 20531, 28548, 25342, 6510, 20657, 15608, 5534, 22145,
|
||||
30150, 25222, 12128, 389, 21860, 9631, 4536, 4704, 3677, 27282, 26668,
|
||||
20784, 15684, 12847, 1307, 10586, 24355, 27553, 10952, 8886, 25029,
|
||||
29278, 29964, 17943, 1006, 5895, 11466, 16679, 17500, 5414, 3420, 17644,
|
||||
5165, 6255, 2807, 30577, 26277, 14032, 2425, 13945, 27988, 17437, 28204,
|
||||
11853, 12265, 8097, 24919, 10703, 18081, 19121, 23364, 14035, 2382,
|
||||
1722, 21617, 11863, 27682, 8538, 26401,
|
||||
1487, 14570, 14213, 18315, 30244, 14611, 25421, 13954, 29802,
|
||||
29118, 5788, 7547, 9710, 21645, 17858, 20672, 2295, 21286, 7217, 30405,
|
||||
5090, 22674, 5747, 5809, 13789, 16385, 23732, 12258, 10944, 14669, 2043,
|
||||
1453, 13510, 12422, 24073, 3025, 28094, 2770, 9198, 27411, 24736, 28958,
|
||||
23508, 27897, 17838, 10690, 5375, 29469, 22458, 9466, 28541, 16308,
|
||||
20491, 10320, 9836, 673, 26630, 20819, 25687, 19263, 16620, 28683,
|
||||
30268, 1113, 26632, 18450, 17555, 20121, 18083, 12796, 26659, 9788,
|
||||
10448, 2828, 29753, 26653, 13636, 6270, 10398, 16224, 1481, 1153, 26387,
|
||||
17835, 19289, 2683, 1937, 16760, 14372, 12632, 15716, 12423, 24202,
|
||||
14543, 10763, 27059, 437, 18647, 17133, 27774,
|
||||
2039, 3931, 7737, 20470, 11068, 26238, 28463, 22610, 28349, 23819,
|
||||
22780, 4101, 13218, 12878, 25048, 25163, 11032, 10129, 2571, 9319,
|
||||
11708, 6704, 19105, 11593, 24863, 26090, 15235, 18038, 22056, 19624,
|
||||
12066, 9798, 16508, 22376, 15776, 10595, 28391, 18898, 11645, 16655,
|
||||
19391, 11364, 28198, 4348, 6653, 11962, 22652, 18750, 22125, 21504,
|
||||
23718, 25662, 6768, 24234, 29605, 8280, 5246, 23064, 1360, 21538, 4374,
|
||||
8186, 7540, 24091, 3017, 23007, 12000, 11289, 8698, 22118, 5505, 18535,
|
||||
29647, 15878, 4416, 8598, 13062, 8878, 9674, 5066, 17770, 24888, 20643,
|
||||
1345, 22570, 1363, 3710, 18429, 11731, 14885, 12983, 18600, 26334,
|
||||
27101, 17858, 22221, 2471, 911, 12033, 2824,
|
||||
6354, 984, 28507, 3521, 17963, 6558, 11166, 24004, 24367, 8572,
|
||||
19198, 6937, 15220, 13122, 3540, 589, 17503, 14073, 14954, 26020, 12974,
|
||||
20684, 19844, 17852, 1097, 10831, 23848, 7013, 15683, 15954, 22290,
|
||||
30257, 15807, 22775, 13607, 9428, 30055, 11607, 30426, 2579, 340, 29747,
|
||||
25213, 28551, 5705, 15704, 10625, 16932, 3215, 16716, 6698, 21470,
|
||||
29839, 511, 23506, 4338, 30506, 18038, 20430, 20586, 18225, 7721, 15812,
|
||||
3140, 22149, 4949, 8125, 9897, 6323, 20612, 2012, 23744, 9414, 16497,
|
||||
5557, 5225, 8518, 30549, 21805, 5692, 25222, 16326, 22995, 27432, 16385,
|
||||
23506, 9911, 23131, 3880, 30647, 13222, 10416, 5619, 2078, 9411, 12398,
|
||||
22772, 7328, 17932, 19965,
|
||||
-1
|
||||
};
|
||||
|
||||
|
||||
int main()
|
||||
{
|
||||
long iter;
|
||||
|
||||
flint_rand_t state;
|
||||
flint_randinit(state);
|
||||
|
||||
printf("cos_minpoly....");
|
||||
fflush(stdout);
|
||||
|
||||
/* compare with table */
|
||||
{
|
||||
long n;
|
||||
fmpz_poly_t p;
|
||||
fmpz_poly_init(p);
|
||||
|
||||
for (n = 0; testdata[n] != -1; n++)
|
||||
{
|
||||
mp_limb_t y;
|
||||
elefun_cos_minpoly(p, n);
|
||||
y = fmpz_poly_evaluate_mod(p, 1337, 31337);
|
||||
|
||||
if (y != testdata[n])
|
||||
{
|
||||
printf("FAIL: n = %ld\n", n);
|
||||
printf("y = %lu\n", y);
|
||||
printf("\n");
|
||||
abort();
|
||||
}
|
||||
}
|
||||
|
||||
fmpz_poly_clear(p);
|
||||
}
|
||||
|
||||
/* check that Phi_n(cos(2pi/n)) = 0 */
|
||||
for (iter = 0; iter < 100; iter++)
|
||||
{
|
||||
long n, prec;
|
||||
fmpz_poly_t phi;
|
||||
fmprb_poly_t fphi;
|
||||
fmprb_t x, y;
|
||||
|
||||
n = 1 + n_randint(state, 2000);
|
||||
prec = 2 + n_randint(state, 1000);
|
||||
|
||||
fmpz_poly_init(phi);
|
||||
fmprb_poly_init(fphi);
|
||||
fmprb_init(x);
|
||||
fmprb_init(y);
|
||||
|
||||
elefun_cos_minpoly(phi, n);
|
||||
fmprb_poly_set_fmpz_poly(fphi, phi, FMPR_PREC_EXACT);
|
||||
|
||||
fmprb_const_pi(x, prec);
|
||||
fmprb_mul_ui(x, x, 2, prec);
|
||||
fmprb_div_ui(x, x, n, prec);
|
||||
fmprb_cos(x, x, prec);
|
||||
|
||||
fmprb_poly_evaluate(y, fphi, x, prec);
|
||||
|
||||
if (!fmprb_contains_zero(y))
|
||||
{
|
||||
printf("FAIL!\n");
|
||||
printf("n = %ld, prec = %ld\n", n, prec);
|
||||
printf("phi = "); fmpz_poly_print(phi); printf("\n\n");
|
||||
printf("x = "); fmprb_printd(x, prec / 3.33); printf("\n\n");
|
||||
printf("y = "); fmprb_printd(y, prec / 3.33); printf("\n\n");
|
||||
abort();
|
||||
}
|
||||
|
||||
fmpz_poly_clear(phi);
|
||||
fmprb_poly_clear(fphi);
|
||||
fmprb_clear(x);
|
||||
fmprb_clear(y);
|
||||
}
|
||||
|
||||
_fmpz_cleanup();
|
||||
printf("PASS\n");
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
|
@ -25,7 +25,7 @@
|
|||
|
||||
#include "fmprb.h"
|
||||
#include "fmprb_poly.h"
|
||||
#include "arith.h"
|
||||
#include "elefun.h"
|
||||
|
||||
void
|
||||
_fmprb_cos_pi_fmpq_algebraic(fmprb_t c, ulong p, ulong q, long prec)
|
||||
|
|
@ -98,9 +98,9 @@ _fmprb_cos_pi_fmpq_algebraic(fmprb_t c, ulong p, ulong q, long prec)
|
|||
fmprb_poly_init(fpoly);
|
||||
|
||||
if (p % 2 == 0)
|
||||
arith_cos_minpoly(poly, q);
|
||||
elefun_cos_minpoly(poly, q);
|
||||
else
|
||||
arith_cos_minpoly(poly, 2 * q);
|
||||
elefun_cos_minpoly(poly, 2 * q);
|
||||
|
||||
eval_extra_prec = fmpz_poly_max_bits(poly);
|
||||
eval_extra_prec = FLINT_ABS(eval_extra_prec);
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue