diff --git a/dev/check_examples.sh b/dev/check_examples.sh index 2f77cd693b..f34c9aeb48 100755 --- a/dev/check_examples.sh +++ b/dev/check_examples.sh @@ -239,6 +239,9 @@ then echo "PASS" exit 0 +elif test "$1" = "minimal_irreducibles"; +then + echo "minimal_irreducibles....SKIPPED" elif test "$1" = "multi_crt"; then echo "multi_crt....SKIPPED" diff --git a/doc/source/fq_nmod.rst b/doc/source/fq_nmod.rst index 0e34f03f23..25f012e45d 100644 --- a/doc/source/fq_nmod.rst +++ b/doc/source/fq_nmod.rst @@ -34,7 +34,7 @@ Context Management Initialises the context for prime `p` and extension degree `d`, with name ``var`` for the generator. By default, it will try - use a Conway polynomial; if one is not available, a random + use a Conway polynomial; if one is not available, a minimal weight irreducible polynomial will be used. Assumes that `p` is a prime. @@ -42,6 +42,17 @@ Context Management Assumes that the string ``var`` is a null-terminated string of length at least one. +.. function:: fq_nmod_ctx_init_minimal_weight_ui(fq_nmod_ctx_t ctx, ulong p, slong d, const char * var) + + Initialises the context for prime `p` and extension degree `d`, + with name ``var`` for the generator, choosing a modulus polynomial + with minimal number of nonzero terms for efficient arithmetic. + + Assumes that `p` is a prime. + + Assumes that the string ``var`` is a null-terminated string + of length at least one. + .. function:: int _fq_nmod_ctx_init_conway_ui(fq_nmod_ctx_t ctx, ulong p, slong d, const char * var) Attempts to initialise the context for prime `p` and extension diff --git a/doc/source/nmod_poly.rst b/doc/source/nmod_poly.rst index 5403788ec5..3adcc74198 100644 --- a/doc/source/nmod_poly.rst +++ b/doc/source/nmod_poly.rst @@ -222,14 +222,47 @@ Randomization Generates a random polynomial with length up to ``len``. -.. function:: void nmod_poly_randtest_irreducible(nmod_poly_t poly, flint_rand_t state, slong len) - - Generates a random irreducible polynomial with length up to ``len``. - .. function:: void nmod_poly_randtest_monic(nmod_poly_t poly, flint_rand_t state, slong len) Generates a random monic polynomial with length ``len``. +.. function:: void nmod_poly_randtest_trinomial(nmod_poly_t poly, flint_rand_t state, slong len) + + Generates a random monic trinomial of length ``len``. + +.. function:: void nmod_poly_randtest_pentomial(nmod_poly_t poly, flint_rand_t state, slong len) + + Generates a random monic pentomial of length ``len``. + +Construction of irreducible polynomials +-------------------------------------------------------------------------------- + +The following functions assume a prime modulus. + +.. function:: void nmod_poly_minimal_irreducible(nmod_poly_t res, ulong n) + + Generates a monic irreducible polynomial of degree ``n`` with minimal + weight (minimal number of nonzero terms). We generate a binomial + if possible, otherwise a trinomial, etc. + It is conjectured that one never needs more than a pentanomial + modulo `p = 2` and a tetranomial modulo `p > 2`. + + More specifically, this function returns the first among all minimal-weight + polynomials in the following ordering. + Firstly, for trinomials, `x^n + a x^k + b` comes before + its monic reversal `x^n + a' x^{n-k} + b'` if `k < n - k`. + Secondly, writing + `f = x^n + a_1 x^{k_1} + a_2 x^{k_2} + \ldots + a_t x^{k_t}` + with `n > k_1 > k_2 \ldots > k_t`, we order tuples + `(a_1, \ldots, a_{t}, k_1, \ldots, k_t)` lexicographically. + We thus favor polynomials with smaller coefficients + (all 1 if possible), and secondly with smaller degrees for the + middle terms. + +.. function:: void nmod_poly_randtest_irreducible(nmod_poly_t poly, flint_rand_t state, slong len) + + Generates a random irreducible polynomial with length up to ``len``. + .. function:: void nmod_poly_randtest_monic_irreducible(nmod_poly_t poly, flint_rand_t state, slong len) Generates a random monic irreducible polynomial with length ``len``. @@ -239,11 +272,6 @@ Randomization Generates a random monic irreducible primitive polynomial with length ``len``. - -.. function:: void nmod_poly_randtest_trinomial(nmod_poly_t poly, flint_rand_t state, slong len) - - Generates a random monic trinomial of length ``len``. - .. function:: int nmod_poly_randtest_trinomial_irreducible(nmod_poly_t poly, flint_rand_t state, slong len, slong max_attempts) Attempts to set ``poly`` to a monic irreducible trinomial of @@ -253,10 +281,6 @@ Randomization trinomials until an irreducible one is found. Returns `1` if one is found and `0` otherwise. -.. function:: void nmod_poly_randtest_pentomial(nmod_poly_t poly, flint_rand_t state, slong len) - - Generates a random monic pentomial of length ``len``. - .. function:: int nmod_poly_randtest_pentomial_irreducible(nmod_poly_t poly, flint_rand_t state, slong len, slong max_attempts) Attempts to set ``poly`` to a monic irreducible pentomial of diff --git a/examples/minimal_irreducibles.c b/examples/minimal_irreducibles.c new file mode 100644 index 0000000000..2fa01b46d7 --- /dev/null +++ b/examples/minimal_irreducibles.c @@ -0,0 +1,103 @@ +/* + Generate irreducible polynomial of minimal weight over GF(p) + for degrees nmin <= n <= nmax. + + This file is public domain. Author: Fredrik Johansson. +*/ + +#include +#include +#include +#include "flint/thread_support.h" +#include "flint/ulong_extras.h" +#include "flint/nmod_poly.h" +#include "flint/nmod_poly_factor.h" +#include "flint/profiler.h" + +typedef struct +{ + nmod_poly_struct * f; + ulong nstart; +} +workinfo_t; + +void +worker(slong i, void * work) +{ + nmod_poly_struct * f = ((workinfo_t *) work)->f; + ulong nstart = ((workinfo_t *) work)->nstart; + nmod_poly_minimal_irreducible(f + i, nstart + i); +} + +int +main(int argc, char * argv[]) +{ + slong i; + int num_threads = 1; + ulong p, n, na, nb, nmin, nmax; + slong batch_size; + workinfo_t work; + + if (argc < 4) + { + flint_printf("usage: minimal_irreducibles [-threads t] p nmin nmax\n"); + return 1; + } + + p = 0; + nmin = 0; + nmax = 0; + + for (i = 1; i < argc; i++) + { + if (!strcmp(argv[i], "-threads")) + { + num_threads = atoi(argv[i+1]); + flint_set_num_threads(num_threads); + i++; + } + else if (i == argc - 3) + p = atol(argv[i]); + else if (i == argc - 2) + nmin = atol(argv[i]); + else if (i == argc - 1) + nmax = atol(argv[i]); + } + + flint_printf("p = %wu, nmin = %wu, nmax = %wu\n", p, nmin, nmax); + + if (!n_is_prime(p)) + flint_throw(FLINT_ERROR, "p = %wu is not prime\n", p); + + if (num_threads == 1) + batch_size = 1; + else + batch_size = FLINT_MIN(10 * num_threads, nmax - nmin + 1); + work.f = flint_malloc(sizeof(nmod_poly_struct) * batch_size); + for (i = 0; i < batch_size; i++) + nmod_poly_init(work.f + i, p); + + TIMEIT_ONCE_START + na = nmin; + while (1) + { + work.nstart = na; + nb = FLINT_MIN(na + batch_size - 1, nmax); + flint_parallel_do(worker, &work, nb - na + 1, 0, 0); + for (n = na; n <= nb; n++) + flint_printf("%{nmod_poly}\n", work.f + n - na); + na = nb + 1; + if (na > nmax) + break; + } + TIMEIT_ONCE_STOP + SHOW_MEMORY_USAGE + + for (i = 0; i < batch_size; i++) + nmod_poly_clear(work.f + i); + flint_free(work.f); + + flint_cleanup_master(); + return 0; +} + diff --git a/src/fq_nmod.h b/src/fq_nmod.h index 8335be0730..26e68917b2 100644 --- a/src/fq_nmod.h +++ b/src/fq_nmod.h @@ -34,6 +34,7 @@ extern "C" { void fq_nmod_ctx_init_ui(fq_nmod_ctx_t ctx, ulong prime, slong deg, const char * var); int _fq_nmod_ctx_init_conway_ui(fq_nmod_ctx_t ctx, ulong prime, slong deg, const char * var); void fq_nmod_ctx_init_conway_ui(fq_nmod_ctx_t ctx, ulong prime, slong deg, const char * var); +void fq_nmod_ctx_init_minimal_weight_ui(fq_nmod_ctx_t ctx, ulong prime, slong deg, const char * var); void fq_nmod_ctx_init_modulus(fq_nmod_ctx_t ctx, const nmod_poly_t modulus, const char * var); void fq_nmod_ctx_init_randtest(fq_nmod_ctx_t ctx, flint_rand_t state, int type); diff --git a/src/fq_nmod/ctx_init.c b/src/fq_nmod/ctx_init.c index 06cccbe9b1..d963918283 100644 --- a/src/fq_nmod/ctx_init.c +++ b/src/fq_nmod/ctx_init.c @@ -51,29 +51,27 @@ void fq_nmod_ctx_init_conway_ui(fq_nmod_ctx_t ctx, ulong p, slong d, const char ctx->is_conway = 1; } +void fq_nmod_ctx_init_minimal_weight_ui(fq_nmod_ctx_t ctx, ulong p, slong d, const char *var) +{ + nmod_poly_t poly; + + ctx->is_conway = 0; + + nmod_poly_init2(poly, p, d + 1); + nmod_poly_minimal_irreducible(poly, d); + fq_nmod_ctx_init_modulus(ctx, poly, var); + nmod_poly_clear(poly); +} + void fq_nmod_ctx_init_ui(fq_nmod_ctx_t ctx, ulong p, slong d, const char *var) { if (_fq_nmod_ctx_init_conway_ui(ctx, p, d, var)) { ctx->is_conway = 1; - return; } else { - flint_rand_t state; - nmod_poly_t poly; - - ctx->is_conway = 0; - - flint_rand_init(state); - - nmod_poly_init2(poly, p, d + 1); - nmod_poly_randtest_sparse_irreducible(poly, state, d + 1); - - fq_nmod_ctx_init_modulus(ctx, poly, var); - - nmod_poly_clear(poly); - flint_rand_clear(state); + fq_nmod_ctx_init_minimal_weight_ui(ctx, p, d, var); } } diff --git a/src/nmod_poly.h b/src/nmod_poly.h index 986a3f9b9f..ca86f6f92a 100644 --- a/src/nmod_poly.h +++ b/src/nmod_poly.h @@ -223,6 +223,10 @@ void nmod_poly_randtest_not_zero(nmod_poly_t poly, flint_rand_t state, slong len } while (nmod_poly_is_zero(poly)); } +/* Construction of irreducibles *********************************************/ + +void nmod_poly_minimal_irreducible(nmod_poly_t res, ulong n); + void nmod_poly_randtest_irreducible(nmod_poly_t poly, flint_rand_t state, slong len); void nmod_poly_randtest_monic(nmod_poly_t poly, flint_rand_t state, slong len); diff --git a/src/nmod_poly/minimal_irreducible.c b/src/nmod_poly/minimal_irreducible.c new file mode 100644 index 0000000000..834f8536aa --- /dev/null +++ b/src/nmod_poly/minimal_irreducible.c @@ -0,0 +1,926 @@ +/* + Copyright (C) 2025 Fredrik Johansson + + This file is part of FLINT. + + FLINT is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. See . +*/ + +#include +#include "ulong_extras.h" +#include "nmod.h" +#include "nmod_vec.h" +#include "nmod_poly.h" +#include "nmod_poly_factor.h" + +static ulong +n_multiplicative_order(ulong x, ulong p, ulong pinv, n_factor_t * p1fac) +{ + ulong m, q, mm; + slong i; + m = p - 1; + + for (i = 0; i < p1fac->num; i++) + { + q = p1fac->p[i]; + + while (m % q == 0) + { + mm = m / q; + if (n_powmod2_preinv(x, mm, p, pinv) != 1) + break; + m = mm; + } + } + + return m; +} + +int +nmod_poly_irreducible_binomial(nmod_poly_t res, ulong n) +{ + ulong q = res->mod.n; + ulong e, a, qinv; + n_factor_t nfac, q1fac; + slong i; + + if (n == 0) + return 0; + + if (n == 1) + { + nmod_poly_one(res); + nmod_poly_set_coeff_ui(res, 1, 1); + return 1; + } + + /* There is an irreducible binomial x^n + a mod q iff + p | q - 1 for each prime p | n and + if 4 | n then 4 | q - 1. */ + if (n % 4 == 0) + if (q % 4 != 1) + return 0; + n_factor_init(&nfac); + n_factor(&nfac, n, 1); + for (i = 0; i < nfac.num; i++) + if ((q - 1) % nfac.p[i] != 0) + return 0; + + nmod_poly_zero(res); + nmod_poly_set_coeff_ui(res, n, 1); + + /* Let e = ord_q(-a). Then we have irreducibility iff + (p divides e but p does not divide (q-1)/e) + for each prime factor p of n. */ + n_factor_init(&q1fac); + n_factor(&q1fac, q - 1, 1); + qinv = n_preinvert_limb(q); + + for (a = 1; a < q; a++) + { + e = n_multiplicative_order(q - a, q, qinv, &q1fac); + + for (i = 0; i < nfac.num; i++) + if (!(e % nfac.p[i] == 0 && ((q - 1) / e) % nfac.p[i] != 0)) + goto next; + + nmod_poly_set_coeff_ui(res, 0, a); + return 1; + + next: + continue; + } + + flint_throw(FLINT_ERROR, "Failed to construct irreducible binomial (n = %wu, p = %wu)", n, q); +} + + + +/* Todo: refine and extend this table */ +static const unsigned int sieve_limit_tab[][10] = { + /* p 1 2 3 4 5 6 7 8 */ + { 2, 1, 1, 1, 1, 24, 24, 256, 256, 0 }, + { 3, 1, 1, 1, 1, 128, 512, 512, 3072, 0 }, + { 5, 1, 1, 1, 48, 512, 1536, 3072, 0 }, + { 7, 1, 1, 1, 64, 1536, 3072, 0 }, + { 11, 1, 1, 64, 768, 2048, 0 }, + { 13, 1, 1, 48, 768, 0 }, + { 17, 1, 1, 384, 768, 0 }, + { 19, 1, 16, 96, 1536, 0 }, + { 23, 1, 16, 256, 2048, 0 }, + { 29, 1, 12, 512, 0 }, + { 31, 1, 12, 384, 0 }, + { 37, 1, 12, 384, 0 }, + { 41, 1, 24, 1024, 0 }, + { 43, 1, 32, 512, 0 }, + { 47, 1, 32, 768, 0 }, + { 53, 1, 48, 1024, 0 }, + { 59, 1, 32, 1024, 0 }, + { 61, 1, 64, 1024, 0 }, + { 67, 1, 64, 1536, 0 }, + { 71, 1, 64, 1536, 0 }, + { 73, 1, 64, 1536, 0 }, + { 79, 1, 64, 1536, 0 }, + { 83, 1, 64, 2048, 0 }, + { 89, 1, 64, 2048, 0 }, + { 97, 4, 128, 2048, 0 }, + { 101, 8, 128, 3072, 0 }, + { 103, 4, 96, 3072, 0 }, + { 107, 4, 96, 0 }, + { 127, 4, 128, 0 }, + { 179, 4, 192, 0 }, + { 269, 8, 256, 0 }, + { 283, 8, 384, 0 }, + { 431, 8, 512, 0 }, + { 487, 8, 768, 0 }, + { 653, 12, 768, 0 }, + { 727, 16, 1024, 0 }, + { 953, 16, 1536, 0 }, + { 1009, 24, 2048, 0 }, + { 1549, 32, 0 }, + { 2333, 48, 0 }, + { 3877, 64, 0 }, + { 7789, 96, 0 }, + { 14107, 128, 0 }, + { 28069, 192, 0 }, + { 35831, 256, 0 }, + { 82457, 384, 0 }, + { 121949, 768, 0 }, + { 180247, 1024, 0 }, + { 0 }, +}; + +/* Return b such that we want to do trial division by all irreducible + factors up to degree b. */ +static ulong trinomial_sieve_limit(ulong p, ulong n) +{ + slong i; + ulong ptab, b, limit; + + for (i = 0; ; i++) + { + ptab = sieve_limit_tab[i][0]; + if (ptab == 0) + return 0; + + if (ptab >= p) + { + limit = 0; + for (b = 1; sieve_limit_tab[i][b] != 0; b++) + { + if (n >= sieve_limit_tab[i][b]) + limit = b; + } + + return limit; + } + } +} + + +/* Return the product of all irreducible polynomial of degree r. */ +/* Well-known fact: x^(p^r) - x contains all factors of degree dividing r, + so we just need to divide out the smaller products. */ +void +nmod_poly_product_all_irreducibles_deg(nmod_poly_t res, ulong r) +{ + ulong k, d, step, p = res->mod.n; + + FLINT_ASSERT(r >= 1 && r <= 8); + + nmod_poly_zero(res); + + if (r == 1) + { + nmod_poly_set_coeff_ui(res, 1, 1); + nmod_poly_set_coeff_ui(res, p, 1); + return; + } + + if (r == 4 || r == 6) + { + d = n_pow(p, r) - n_pow(p, 2); + step = p * p - 1; + } + else if (r == 8) + { + d = n_pow(p, r) - n_pow(p, 4); + step = n_pow(p, 4) - 1; + } + else + { + d = n_pow(p, r) - p; + step = p - 1; + } + + for (k = 0; k <= d; k += step) + nmod_poly_set_coeff_ui(res, k, 1); + + if (r == 6) + { + nmod_poly_t t; + nmod_poly_init(t, p); + nmod_poly_product_all_irreducibles_deg(t, 3); + nmod_poly_divexact(res, res, t); + nmod_poly_clear(t); + } +} + +/* Hack: nmod_poly_gcd is not optimised for unbalanced GCD with the smaller + operand sparse, so we deal with this here. */ +void +_nmod_poly_inplace_rem_sparse_monic(nn_ptr R, + slong lenA, nn_srcptr Bcoeffs, const slong * Bexps, slong nzB, slong lenB, nmod_t mod) +{ + slong i, j, k; + slong n = lenB - 1; + ulong c; + + for (i = lenA - 1; i >= n; i--) + { + c = R[i]; + + for (k = nzB - 2; k >= 0; k--) + { + j = Bexps[k]; + R[j + i - n] = nmod_sub(R[j + i - n], nmod_mul(c, Bcoeffs[k], mod), mod); + } + } +} + +void +nmod_poly_gcd_with_sparse(nmod_poly_t res, const nmod_poly_t A, const nmod_poly_t B, + nn_srcptr Bcoeffs, const slong * Bexps, slong nzB) +{ + if (A->length <= B->length) + { + nmod_poly_gcd(res, A, B); + } + else + { + nn_ptr R; + slong Rlen = B->length - 1; + nmod_poly_t T; + TMP_INIT; + TMP_START; + R = TMP_ALLOC(sizeof(ulong) * A->length); + _nmod_vec_set(R, A->coeffs, A->length); + _nmod_poly_inplace_rem_sparse_monic(R, A->length, Bcoeffs, Bexps, nzB, B->length, res->mod); + NMOD_VEC_NORM(R, Rlen); + T->coeffs = R; + T->length = Rlen; + T->alloc = Rlen; + nmod_poly_gcd(res, T, B); + TMP_END; + } +} + +/* Discriminant of x^n + ax^k + b mod p > 2 (Swan's formula). */ +ulong +_nmod_poly_trinomial_discriminant(ulong n, ulong k, ulong a, ulong b, nmod_t mod) +{ + ulong d, n1, k1; + ulong D1, D2, D3; + + d = n_gcd(n, k); + n1 = n / d; + k1 = k / d; + + D1 = nmod_pow_ui(b, k - 1, mod); + D2 = nmod_mul(nmod_pow_ui(n, n1, mod), nmod_pow_ui(b, n1 - k1, mod), mod); + D3 = nmod_mul(nmod_pow_ui(n - k, n1 - k1, mod), nmod_pow_ui(k, k1, mod), mod); + D3 = nmod_mul(D3, nmod_pow_ui(a, n1, mod), mod); + if (n1 % 2) + D3 = nmod_neg(D3, mod); + D2 = nmod_sub(D2, D3, mod); + D2 = nmod_pow_ui(D2, d, mod); + D1 = nmod_mul(D1, D2, mod); + if ((n * (n - 1) / 2) % 2) + D1 = nmod_neg(D1, mod); + + return D1; +} + +/* Exceptional n for which no irreducible trinomial exists. These tables are + not needed for correctness, but speed up the search. */ + +/* All exceptions for n <= 10000 */ +static const uint16_t no_trinomial_mod3[] = { + 49, 57, 65, 68, 75, 98, 105, 123, 129, 130, 132, 149, 161, 175, 189, 197, + 207, 212, 213, 221, 223, 231, 233, 264, 267, 276, 281, 292, 297, 298, 303, 309, + 311, 319, 332, 343, 391, 394, 397, 401, 404, 405, 411, 415, 426, 435, 436, 437, + 439, 441, 442, 453, 459, 463, 489, 497, 513, 521, 524, 528, 534, 540, 543, 545, + 548, 551, 552, 556, 559, 561, 571, 572, 591, 597, 612, 613, 619, 622, 627, 628, + 633, 636, 641, 653, 657, 664, 671, 683, 686, 692, 699, 703, 708, 715, 723, 725, + 732, 737, 747, 748, 753, 761, 783, 794, 795, 799, 802, 804, 809, 822, 823, 825, + 836, 847, 851, 852, 855, 874, 878, 879, 881, 891, 893, 903, 905, 908, 909, 918, + 933, 951, 953, 957, 961, 965, 975, 979, 981, 989, 994, 996, 1003, 1005, 1011, 1029, + 1039, 1041, 1047, 1051, 1056, 1059, 1065, 1068, 1073, 1075, 1080, 1086, 1097, 1099, 1102, 1104, + 1109, 1111, 1125, 1127, 1129, 1135, 1140, 1142, 1147, 1153, 1164, 1193, 1195, 1203, 1221, 1224, + 1225, 1226, 1227, 1233, 1238, 1239, 1241, 1252, 1253, 1261, 1263, 1267, 1269, 1272, 1277, 1281, + 1297, 1308, 1316, 1342, 1345, 1356, 1357, 1363, 1365, 1366, 1385, 1389, 1393, 1411, 1416, 1419, + 1431, 1435, 1445, 1450, 1452, 1464, 1467, 1471, 1481, 1484, 1491, 1500, 1505, 1507, 1523, 1524, + 1533, 1555, 1565, 1566, 1572, 1581, 1587, 1589, 1603, 1605, 1608, 1611, 1617, 1623, 1628, 1639, + 1644, 1646, 1647, 1649, 1673, 1677, 1687, 1689, 1694, 1695, 1697, 1704, 1709, 1713, 1723, 1725, + 1731, 1733, 1747, 1748, 1749, 1762, 1785, 1786, 1788, 1805, 1810, 1829, 1831, 1836, 1851, 1869, + 1873, 1876, 1877, 1883, 1887, 1889, 1892, 1902, 1905, 1906, 1914, 1917, 1922, 1925, 1930, 1937, + 1957, 1958, 1967, 1971, 1975, 1978, 1992, 2004, 2005, 2007, 2010, 2019, 2023, 2031, 2047, 2057, + 2061, 2071, 2078, 2103, 2108, 2112, 2116, 2118, 2124, 2130, 2136, 2139, 2143, 2146, 2150, 2153, + 2155, 2160, 2172, 2181, 2189, 2198, 2203, 2204, 2205, 2208, 2213, 2217, 2218, 2220, 2222, 2229, + 2241, 2254, 2257, 2258, 2259, 2265, 2270, 2273, 2280, 2283, 2284, 2297, 2299, 2308, 2311, 2316, + 2323, 2328, 2331, 2333, 2339, 2348, 2357, 2367, 2369, 2375, 2386, 2390, 2395, 2405, 2413, 2417, + 2428, 2429, 2433, 2441, 2448, 2451, 2454, 2455, 2457, 2460, 2469, 2476, 2478, 2479, 2487, 2489, + 2501, 2511, 2517, 2526, 2527, 2534, 2538, 2541, 2544, 2554, 2555, 2557, 2559, 2561, 2571, 2591, + 2595, 2611, 2616, 2633, 2636, 2643, 2644, 2649, 2657, 2668, 2671, 2681, 2684, 2690, 2700, 2708, + 2711, 2712, 2715, 2733, 2745, 2755, 2769, 2770, 2773, 2779, 2781, 2786, 2787, 2791, 2801, 2803, + 2809, 2815, 2825, 2827, 2831, 2832, 2838, 2843, 2844, 2862, 2863, 2870, 2890, 2904, 2908, 2909, + 2921, 2923, 2924, 2928, 2935, 2942, 2955, 2959, 2961, 2962, 2965, 2993, 3000, 3003, 3010, 3014, + 3017, 3027, 3029, 3033, 3043, 3048, 3049, 3051, 3053, 3060, 3063, 3067, 3069, 3076, 3077, 3087, + 3091, 3093, 3116, 3117, 3119, 3123, 3124, 3125, 3127, 3132, 3140, 3144, 3151, 3165, 3167, 3171, + 3178, 3183, 3185, 3204, 3206, 3207, 3212, 3216, 3219, 3222, 3228, 3244, 3246, 3249, 3252, 3261, + 3269, 3281, 3288, 3291, 3292, 3297, 3308, 3329, 3332, 3346, 3348, 3357, 3367, 3369, 3387, 3396, + 3401, 3408, 3413, 3418, 3423, 3428, 3429, 3435, 3439, 3446, 3449, 3459, 3466, 3471, 3481, 3487, + 3507, 3508, 3540, 3556, 3567, 3572, 3576, 3585, 3593, 3596, 3597, 3603, 3605, 3609, 3610, 3615, + 3619, 3621, 3643, 3662, 3663, 3668, 3672, 3684, 3687, 3689, 3691, 3703, 3711, 3716, 3717, 3724, + 3737, 3746, 3749, 3751, 3754, 3756, 3764, 3765, 3766, 3774, 3777, 3778, 3780, 3787, 3791, 3797, + 3804, 3810, 3825, 3828, 3837, 3844, 3853, 3855, 3859, 3860, 3861, 3871, 3875, 3876, 3885, 3891, + 3909, 3914, 3919, 3929, 3931, 3932, 3934, 3941, 3945, 3950, 3953, 3955, 3963, 3972, 3977, 3981, + 3984, 3999, 4004, 4008, 4010, 4017, 4020, 4023, 4025, 4027, 4029, 4039, 4041, 4062, 4065, 4073, + 4119, 4122, 4123, 4131, 4133, 4135, 4137, 4140, 4142, 4143, 4156, 4157, 4161, 4171, 4197, 4216, + 4217, 4219, 4224, 4227, 4232, 4236, 4244, 4245, 4248, 4257, 4260, 4261, 4272, 4278, 4286, 4289, + 4293, 4305, 4308, 4313, 4315, 4320, 4327, 4344, 4362, 4363, 4380, 4385, 4389, 4396, 4416, 4426, + 4433, 4436, 4440, 4444, 4445, 4447, 4452, 4457, 4459, 4469, 4471, 4473, 4476, 4479, 4487, 4493, + 4495, 4499, 4503, 4508, 4519, 4523, 4530, 4531, 4540, 4545, 4546, 4547, 4548, 4560, 4565, 4569, + 4571, 4581, 4591, 4594, 4597, 4601, 4613, 4622, 4627, 4632, 4635, 4646, 4647, 4651, 4652, 4653, + 4656, 4663, 4665, 4666, 4668, 4678, 4684, 4697, 4707, 4711, 4714, 4721, 4738, 4747, 4750, 4753, + 4771, 4772, 4773, 4781, 4791, 4796, 4817, 4833, 4834, 4836, 4841, 4843, 4855, 4858, 4860, 4876, + 4882, 4896, 4899, 4901, 4902, 4905, 4907, 4908, 4910, 4920, 4923, 4935, 4938, 4949, 4956, 4971, + 4975, 4985, 4987, 4989, 4993, 4997, 5002, 5009, 5011, 5015, 5017, 5023, 5041, 5045, 5052, 5054, + 5055, 5059, 5061, 5068, 5069, 5076, 5079, 5083, 5084, 5087, 5088, 5093, 5097, 5108, 5122, 5132, + 5151, 5155, 5165, 5171, 5182, 5193, 5199, 5211, 5212, 5219, 5222, 5228, 5232, 5257, 5261, 5268, + 5286, 5297, 5301, 5303, 5313, 5314, 5321, 5333, 5342, 5345, 5349, 5357, 5361, 5362, 5369, 5372, + 5373, 5391, 5393, 5400, 5416, 5417, 5422, 5424, 5430, 5435, 5453, 5457, 5465, 5466, 5481, 5487, + 5491, 5492, 5509, 5510, 5513, 5515, 5524, 5525, 5538, 5539, 5541, 5546, 5553, 5561, 5565, 5567, + 5571, 5575, 5582, 5588, 5601, 5602, 5606, 5612, 5618, 5621, 5630, 5635, 5652, 5655, 5659, 5662, + 5664, 5673, 5676, 5679, 5683, 5688, 5697, 5703, 5707, 5713, 5721, 5724, 5726, 5733, 5739, 5755, + 5763, 5764, 5765, 5780, 5787, 5791, 5799, 5801, 5804, 5808, 5818, 5823, 5825, 5828, 5837, 5839, + 5842, 5844, 5846, 5852, 5856, 5870, 5876, 5892, 5895, 5905, 5910, 5918, 5919, 5923, 5925, 5929, + 5930, 5931, 5933, 5935, 5941, 5945, 5955, 5957, 5967, 5969, 5989, 5995, 6000, 6003, 6005, 6029, + 6031, 6034, 6036, 6039, 6052, 6054, 6060, 6066, 6079, 6096, 6099, 6101, 6102, 6106, 6111, 6116, + 6120, 6134, 6137, 6141, 6161, 6163, 6165, 6177, 6180, 6186, 6191, 6207, 6211, 6212, 6223, 6243, + 6250, 6251, 6254, 6255, 6259, 6264, 6265, 6283, 6287, 6288, 6301, 6317, 6329, 6330, 6334, 6356, + 6361, 6363, 6365, 6370, 6372, 6377, 6379, 6393, 6408, 6412, 6414, 6425, 6427, 6432, 6444, 6456, + 6463, 6465, 6473, 6476, 6477, 6483, 6485, 6489, 6492, 6493, 6495, 6504, 6508, 6524, 6531, 6547, + 6548, 6551, 6559, 6567, 6576, 6579, 6585, 6587, 6603, 6609, 6619, 6621, 6629, 6652, 6675, 6677, + 6687, 6689, 6696, 6714, 6716, 6734, 6738, 6739, 6743, 6772, 6775, 6783, 6792, 6795, 6802, 6811, + 6816, 6826, 6847, 6856, 6858, 6867, 6868, 6893, 6898, 6905, 6909, 6919, 6937, 6948, 6949, 6955, + 6957, 6967, 6974, 6999, 7005, 7011, 7025, 7037, 7041, 7067, 7075, 7079, 7080, 7087, 7093, 7099, + 7111, 7113, 7121, 7123, 7141, 7149, 7152, 7161, 7167, 7181, 7183, 7186, 7197, 7203, 7207, 7219, + 7221, 7229, 7230, 7231, 7235, 7238, 7241, 7252, 7263, 7267, 7275, 7286, 7293, 7301, 7305, 7308, + 7313, 7339, 7344, 7345, 7348, 7349, 7359, 7363, 7364, 7368, 7373, 7380, 7382, 7383, 7385, 7395, + 7406, 7407, 7412, 7429, 7433, 7435, 7452, 7467, 7471, 7474, 7483, 7491, 7493, 7495, 7500, 7512, + 7530, 7532, 7533, 7548, 7553, 7556, 7559, 7560, 7567, 7575, 7594, 7597, 7604, 7608, 7613, 7617, + 7620, 7625, 7628, 7629, 7635, 7637, 7639, 7641, 7649, 7653, 7656, 7663, 7669, 7673, 7683, 7687, + 7688, 7692, 7697, 7706, 7707, 7710, 7713, 7717, 7727, 7731, 7732, 7735, 7737, 7742, 7745, 7750, + 7752, 7755, 7757, 7761, 7769, 7775, 7779, 7797, 7803, 7818, 7819, 7827, 7833, 7838, 7853, 7855, + 7857, 7858, 7862, 7867, 7876, 7881, 7882, 7889, 7890, 7895, 7903, 7924, 7926, 7944, 7945, 7948, + 7951, 7959, 7965, 7968, 7971, 7972, 7973, 7980, 7981, 7993, 7997, 7998, 8004, 8009, 8016, 8021, + 8023, 8027, 8033, 8037, 8040, 8050, 8052, 8057, 8059, 8068, 8075, 8078, 8079, 8081, 8084, 8091, + 8095, 8097, 8100, 8103, 8109, 8119, 8124, 8129, 8132, 8143, 8146, 8155, 8156, 8167, 8179, 8199, + 8201, 8203, 8205, 8209, 8215, 8220, 8238, 8239, 8244, 8246, 8251, 8253, 8259, 8261, 8262, 8265, + 8266, 8268, 8270, 8271, 8273, 8276, 8280, 8284, 8289, 8321, 8323, 8333, 8342, 8343, 8356, 8371, + 8388, 8396, 8405, 8417, 8429, 8434, 8436, 8438, 8441, 8448, 8451, 8453, 8468, 8469, 8472, 8475, + 8496, 8513, 8514, 8520, 8522, 8527, 8532, 8533, 8544, 8549, 8553, 8556, 8557, 8571, 8573, 8575, + 8586, 8589, 8597, 8609, 8610, 8616, 8626, 8630, 8640, 8643, 8647, 8649, 8657, 8668, 8673, 8681, + 8685, 8688, 8691, 8693, 8695, 8700, 8709, 8717, 8724, 8726, 8733, 8743, 8757, 8760, 8764, 8767, + 8769, 8772, 8779, 8781, 8789, 8799, 8801, 8803, 8805, 8807, 8823, 8825, 8827, 8832, 8835, 8844, + 8845, 8847, 8849, 8852, 8857, 8865, 8880, 8889, 8890, 8894, 8904, 8911, 8913, 8914, 8918, 8923, + 8925, 8927, 8935, 8937, 8940, 8942, 8943, 8949, 8952, 8958, 8961, 8971, 8974, 8985, 8986, 8990, + 8993, 8997, 8998, 9004, 9006, 9012, 9019, 9021, 9033, 9039, 9041, 9046, 9057, 9060, 9065, 9069, + 9092, 9094, 9096, 9099, 9103, 9105, 9116, 9117, 9120, 9124, 9127, 9137, 9138, 9142, 9145, 9147, + 9148, 9153, 9156, 9163, 9180, 9182, 9187, 9189, 9194, 9213, 9223, 9225, 9235, 9239, 9254, 9261, + 9264, 9270, 9285, 9294, 9302, 9305, 9309, 9312, 9317, 9329, 9330, 9335, 9336, 9340, 9341, 9351, + 9353, 9367, 9372, 9377, 9387, 9389, 9394, 9396, 9399, 9405, 9407, 9411, 9413, 9414, 9420, 9422, + 9433, 9435, 9436, 9444, 9447, 9452, 9471, 9476, 9492, 9497, 9499, 9506, 9507, 9509, 9513, 9524, + 9529, 9540, 9542, 9544, 9545, 9546, 9549, 9559, 9561, 9562, 9569, 9571, 9579, 9582, 9587, 9595, + 9596, 9597, 9605, 9612, 9615, 9617, 9619, 9621, 9637, 9639, 9644, 9653, 9655, 9665, 9666, 9669, + 9672, 9677, 9682, 9684, 9686, 9697, 9701, 9705, 9708, 9710, 9711, 9716, 9720, 9737, 9749, 9759, + 9773, 9779, 9780, 9783, 9792, 9796, 9798, 9802, 9804, 9807, 9809, 9813, 9816, 9819, 9831, 9840, + 9847, 9849, 9855, 9875, 9876, 9877, 9879, 9881, 9893, 9901, 9907, 9911, 9912, 9921, 9924, 9931, + 9933, 9943, 9950, 9951, 9953, 9957, 9970, 9975, 9977, 9978, 9979, 9980, 0 +}; + +/* All exceptions for n <= 5000 */ +static const uint16_t no_trinomial_mod5[] = { + 35, 70, 123, 125, 140, 181, 191, 209, 213, 219, 237, 249, 250, 253, 265, 273, + 280, 285, 307, 319, 345, 351, 362, 375, 382, 407, 413, 415, 418, 421, 429, 433, + 441, 445, 457, 461, 498, 500, 519, 530, 539, 547, 560, 570, 581, 638, 663, 683, + 690, 699, 723, 724, 725, 735, 750, 753, 765, 779, 787, 795, 803, 815, 817, 830, + 833, 835, 859, 866, 885, 890, 915, 917, 927, 957, 987, 993, 1000, 1015, 1037, 1038, + 1043, 1055, 1060, 1087, 1095, 1113, 1120, 1121, 1137, 1140, 1159, 1161, 1173, 1183, 1185, 1203, + 1249, 1257, 1261, 1267, 1275, 1276, 1293, 1299, 1335, 1351, 1357, 1359, 1365, 1380, 1381, 1387, + 1407, 1427, 1443, 1449, 1450, 1470, 1473, 1479, 1500, 1501, 1503, 1515, 1517, 1530, 1537, 1557, + 1558, 1575, 1590, 1606, 1619, 1621, 1630, 1631, 1653, 1660, 1665, 1670, 1671, 1672, 1683, 1687, + 1691, 1718, 1719, 1727, 1732, 1749, 1769, 1770, 1780, 1787, 1827, 1830, 1831, 1833, 1834, 1837, + 1855, 1869, 1873, 1877, 1887, 1899, 1901, 1911, 1914, 1929, 1959, 1967, 1977, 1986, 2000, 2025, + 2030, 2037, 2041, 2076, 2087, 2091, 2097, 2110, 2120, 2141, 2151, 2157, 2163, 2173, 2187, 2190, + 2240, 2242, 2251, 2257, 2261, 2280, 2319, 2327, 2339, 2341, 2345, 2353, 2359, 2367, 2370, 2379, + 2383, 2407, 2421, 2433, 2447, 2459, 2463, 2485, 2491, 2517, 2519, 2522, 2535, 2541, 2550, 2552, + 2575, 2583, 2595, 2601, 2607, 2621, 2631, 2645, 2653, 2670, 2673, 2683, 2702, 2705, 2717, 2730, + 2745, 2755, 2757, 2760, 2762, 2763, 2775, 2809, 2843, 2847, 2867, 2871, 2877, 2879, 2885, 2891, + 2895, 2899, 2900, 2917, 2940, 2967, 2991, 3000, 3002, 3003, 3027, 3030, 3035, 3049, 3051, 3053, + 3060, 3074, 3095, 3103, 3125, 3131, 3150, 3165, 3173, 3180, 3203, 3212, 3213, 3242, 3260, 3275, + 3295, 3320, 3327, 3330, 3340, 3344, 3361, 3366, 3377, 3387, 3395, 3403, 3435, 3439, 3447, 3455, + 3463, 3464, 3469, 3498, 3525, 3527, 3540, 3553, 3555, 3560, 3574, 3645, 3660, 3675, 3687, 3710, + 3715, 3723, 3727, 3746, 3759, 3795, 3802, 3803, 3807, 3823, 3831, 3877, 3885, 3889, 3895, 3901, + 3907, 3911, 3925, 3943, 3945, 3951, 3953, 3969, 3972, 3984, 4000, 4013, 4029, 4041, 4050, 4060, + 4077, 4079, 4095, 4111, 4137, 4147, 4151, 4163, 4165, 4174, 4211, 4220, 4227, 4240, 4263, 4282, + 4293, 4305, 4325, 4346, 4367, 4380, 4391, 4409, 4425, 4443, 4465, 4479, 4480, 4484, 4487, 4489, + 4497, 4503, 4560, 4563, 4573, 4605, 4635, 4647, 4655, 4663, 4682, 4687, 4690, 4691, 4697, 4706, + 4740, 4755, 4767, 4803, 4807, 4825, 4839, 4851, 4879, 4894, 4901, 4905, 4911, 4915, 4918, 4947, + 4951, 4953, 4970, 4981, 4995, 4997, 0 +}; + +/* All exceptions for n <= 5000 */ +static const uint16_t no_trinomial_mod7[] = { + 124, 163, 268, 301, 359, 364, 372, 385, 389, 407, 431, 476, 489, 527, 563, 601, + 628, 649, 728, 829, 835, 853, 862, 952, 995, 1031, 1063, 1077, 1079, 1115, 1259, 1285, + 1343, 1355, 1372, 1384, 1444, 1456, 1467, 1477, 1669, 1703, 1705, 1706, 1735, 1777, 1781, 1803, + 1825, 1828, 1884, 1889, 1904, 1975, 2003, 2044, 2057, 2083, 2093, 2105, 2158, 2189, 2281, 2339, + 2395, 2505, 2507, 2513, 2564, 2668, 2744, 2755, 2767, 2773, 2807, 2812, 2833, 2863, 2867, 2912, + 2985, 3031, 3091, 3093, 3157, 3203, 3263, 3271, 3289, 3332, 3345, 3406, 3463, 3487, 3515, 3563, + 3650, 3656, 3716, 3719, 3763, 3797, 3808, 3836, 3875, 3892, 3907, 3916, 3971, 4045, 4055, 4088, + 4116, 4166, 4204, 4210, 4249, 4321, 4340, 4391, 4435, 4607, 4673, 4681, 4687, 4708, 4748, 4763, + 4780, 4781, 4961, 0 +}; + +/* We skip the mod 2 table. Although the exceptions are quite dense, + the performance impact is barely noticeable (most time will be spent + searching pentanomials). It would be more worthwhile to just include + a table of the irreducible polynomials themselves over GF(2). */ + +/* Exceptions are also known mod 11, 13 and 17, but they are very sparse + and hence not really worth it. */ + +/* +static const uint16_t no_trinomial_mod2[] = { + 8, 13, 16, 19, 24, 26, 27, 32, 37, 38, 40, 43, 45, 48, 50, 51, + 53, 56, 59, 61, 64, 67, 69, 70, 72, 75, 77, 78, 80, 82, 83, 85, + // ... +}; + +// All exceptions for n <= 5000 +static const uint16_t no_trinomial_mod11[] = { + 219, 291, 467, 599, 761, 909, 1137, 1243, 2859, 2873, 2883, 3047, 3531, 3651, 3932, 4139, + 4207, 4247, 4383, 4419, 4508, 4545, 4551, 4564, 4741, 0 +}; + +// All exceptions for n <= 2000 +static const uint16_t no_trinomial_mod13[] = { + 833, 0 +}; + +// All exceptions for n <= 2000 +static const uint16_t no_trinomial_mod17[] = { + 231, 375, 675, 717, 1941, 0 +}; */ + +int +nmod_poly_irreducible_trinomial(nmod_poly_t res, ulong n) +{ + ulong p = res->mod.n; + nmod_t mod = res->mod; + ulong a, b, c, k, l; + ulong a2, b2, e2n, cn, ck, cn2, ck2, D; + slong i; + nn_ptr invctab = NULL; + int want_pruning; + int found = 0; + nmod_poly_t sieve[10], h; + slong num_sieve = 0; + ulong Bcoeffs[3]; + slong Bexps[3]; + ulong sieve_limit; + + if (n <= 1) + return 0; + + if (p == 2 && (n % 8 == 0)) + return 0; + + /* Ciet, Quisquater and Sica, A Short Note on Irreducible Trinomials in Binary Fields */ + if (p == 2 && ((n % 24 == 13) || (n % 24 == 19)) && n_is_prime(n)) + return 0; + + if (p == 3 || p == 5 || p == 7) + { + const uint16_t * tab; + + if (p == 3) tab = no_trinomial_mod3; + else if (p == 5) tab = no_trinomial_mod5; + else tab = no_trinomial_mod7; + + for (i = 0; tab[i] != 0 && tab[i] <= n; i++) + if (tab[i] == n) + return 0; + } + + sieve_limit = trinomial_sieve_limit(p, n); + want_pruning = (sieve_limit > 0); + + if (want_pruning) + { + invctab = _nmod_vec_init(p); + for (c = 1; c < p; c++) + invctab[c] = nmod_inv(c, mod); + + nmod_poly_init_mod(h, mod); + + for (a = 2; a <= sieve_limit && a < n; a++) + { + nmod_poly_init_mod(sieve[num_sieve], mod); + nmod_poly_product_all_irreducibles_deg(sieve[num_sieve], a); + num_sieve++; + } + } + +#define ALREADY_CHECKED(a,b,a2,b2) ((a2) < (a) || ((a2) == (a) && (b2) < (b))) + + /* Iterate over x^n + a x^k + b */ + for (a = 1; a < p; a++) + { + for (b = 1; b < p; b++) + { + /* Mirror symmetry: f(x) is irreducible iff x^n f(1/x) is irreducible */ + for (k = 1; k * 2 <= n; k++) + { + if (p > 2) + { + /* See if we have already checked f(-x). */ + /* n even, k even x^n + ax^k + b -> x^n + ax^k + b */ + /* n even, k odd x^n + ax^k + b -> x^n - ax^k + b */ + /* n odd, k even x^n + ax^k + b -> x^n - ax^k - b */ + /* n odd, k odd x^n + ax^k + b -> x^n + ax^k - b */ + b2 = (n % 2) ? p - b : b; + a2 = ((n ^ k) % 2) ? p - a : a; + if (ALREADY_CHECKED(a, b, a2, b2)) + goto next_trinomial; + } + + /* Swan-Stickelberger */ + if (p == 2) + { + /* Perfect square */ + if (n % 2 == 0 && k % 2 == 0) + goto next_trinomial; + + if (n % 2 == k % 2) + l = n - k; + else + l = k; + + if (n % 2 == 0 && l % 2 == 1 && n != 2 * l && ((n * l) / 2) % 4 <= 1) + goto next_trinomial; + if (n % 2 == 1 && l % 2 == 0 && ((2*n) % l != 0) && (n % 8 == 3 || n % 8 == 5)) + goto next_trinomial; + if (n % 2 == 1 && l % 2 == 0 && ((2*n) % l == 0) && (n % 8 == 1 || n % 8 == 7)) + goto next_trinomial; + } + else + { + D = _nmod_poly_trinomial_discriminant(n, k, a, b, mod); + + if (nmod_pow_ui(D, (p - 1) / 2, mod) == ((n % 2 == 1) ? p - 1 : 1)) + goto next_trinomial; + } + + if (want_pruning && p > 2) + { + for (c = 1; 2 * c < p; c++) + { + ck = nmod_pow_ui(c, k, mod); + cn = nmod_mul(ck, nmod_pow_ui(c, n - k, mod), mod); + + if (c >= 2) + { + /* See if we have already checked f(c*x). */ + e2n = invctab[cn]; + b2 = nmod_mul(b, e2n, mod); + a2 = nmod_mul(nmod_mul(a, ck, mod), e2n, mod); + if (ALREADY_CHECKED(a, b, a2, b2)) + goto next_trinomial; + + /* See if we have already checked f(-c*x). */ + b2 = (n % 2) ? nmod_neg(b2, mod) : b2; + a2 = ((n ^ k) % 2) ? nmod_neg(a2, mod) : a2; + if (ALREADY_CHECKED(a, b, a2, b2)) + goto next_trinomial; + } + + /* Trial division: check if c^n + a c^k + b = 0. */ + if (nmod_add(cn, nmod_mul(a, ck, mod), mod) == p - b) + goto next_trinomial; + + /* Trial division: check if (-c)^n + a (-c)^k + b = 0. */ + cn2 = (n % 2) ? nmod_neg(cn, mod) : cn; + ck2 = (k % 2) ? nmod_neg(ck, mod) : ck; + if (nmod_add(cn2, nmod_mul(a, ck2, mod), mod) == p - b) + goto next_trinomial; + } + } + + nmod_poly_zero(res); + nmod_poly_set_coeff_ui(res, n, 1); + nmod_poly_set_coeff_ui(res, 0, b); + nmod_poly_set_coeff_ui(res, k, a); + + for (i = 0; i < num_sieve; i++) + { + Bcoeffs[0] = b; + Bcoeffs[1] = a; + Bcoeffs[2] = 1; + Bexps[0] = 0; + Bexps[1] = k; + Bexps[2] = n; + + nmod_poly_gcd_with_sparse(h, sieve[i], res, Bcoeffs, Bexps, 3); + if (!nmod_poly_is_one(h)) + goto next_trinomial; + } + + if (nmod_poly_is_irreducible_ddf(res)) + { + found = 1; + goto cleanup; + } + + next_trinomial: + continue; + } + } + } + +cleanup: + if (want_pruning) + { + _nmod_vec_clear(invctab); + for (i = 0; i < num_sieve; i++) + nmod_poly_clear(sieve[i]); + nmod_poly_clear(h); + } + + return found; +} + +int +nmod_poly_irreducible_tetranomial(nmod_poly_t res, ulong n) +{ + ulong p = res->mod.n; + nmod_t mod = res->mod; + ulong a, b, c, d, k, l; + ulong dk, dl, dn, dk2, dl2, dn2, s; + slong i; + int want_pruning; + int found = 0; + nmod_poly_t sieve[10], h; + slong num_sieve = 0; + ulong Bcoeffs[4]; + slong Bexps[4]; + ulong sieve_limit; + + if (n <= 2 || p == 2) + return 0; + + sieve_limit = trinomial_sieve_limit(p, n); + want_pruning = (sieve_limit > 0); + + if (want_pruning) + { + nmod_poly_init_mod(h, mod); + + for (a = 2; a <= sieve_limit && a < n; a++) + { + nmod_poly_init_mod(sieve[num_sieve], mod); + nmod_poly_product_all_irreducibles_deg(sieve[num_sieve], a); + num_sieve++; + } + } + + /* Iterate over all x^n + a x^k + b x^l + c */ + for (a = 1; a < p; a++) + { + for (b = 1; b < p; b++) + { + for (c = 1; c < p; c++) + { + for (k = 2; k < n; k++) + { + for (l = 1; l < k; l++) + { + if (want_pruning) + { + for (d = 1; 2 * d < p; d++) + { + dk = nmod_pow_ui(d, k, mod); + dl = nmod_pow_ui(d, l, mod); + dn = nmod_pow_ui(d, n, mod); + + /* Trial division: check if d^n + a d^k + b d^l + c = 0. */ + s = nmod_add(dn, nmod_mul(a, dk, mod), mod); + s = nmod_add(s, nmod_mul(b, dl, mod), mod); + if (s == p - c) + goto next_tetranomial; + + dn2 = (n % 2) ? nmod_neg(dn, mod) : dn; + dk2 = (k % 2) ? nmod_neg(dk, mod) : dk; + dl2 = (l % 2) ? nmod_neg(dl, mod) : dl; + /* Trial division: check if (-d)^n + a (-d)^k + b (-d)^l + c = 0. */ + s = nmod_add(dn2, nmod_mul(a, dk2, mod), mod); + s = nmod_add(s, nmod_mul(b, dl2, mod), mod); + if (s == p - c) + goto next_tetranomial; + } + } + + nmod_poly_zero(res); + nmod_poly_set_coeff_ui(res, n, 1); + nmod_poly_set_coeff_ui(res, k, a); + nmod_poly_set_coeff_ui(res, l, b); + nmod_poly_set_coeff_ui(res, 0, c); + + for (i = 0; i < num_sieve; i++) + { + Bcoeffs[0] = c; + Bcoeffs[1] = b; + Bcoeffs[2] = a; + Bcoeffs[3] = 1; + Bexps[0] = 0; + Bexps[1] = l; + Bexps[2] = k; + Bexps[3] = n; + + nmod_poly_gcd_with_sparse(h, sieve[i], res, Bcoeffs, Bexps, 4); + if (!nmod_poly_is_one(h)) + goto next_tetranomial; + } + + { + ulong D, D1, D2; + D = nmod_poly_discriminant(res); + D1 = nmod_pow_ui(D, (p - 1) / 2, mod); + D2 = (n % 2 == 1) ? p - 1 : 1; + if (D1 == D2) + goto next_tetranomial; + } + + if (nmod_poly_is_irreducible_ddf(res)) + { + found = 1; + goto cleanup; + } + + next_tetranomial: + continue; + } + } + } + } + } + +cleanup: + if (want_pruning) + { + for (i = 0; i < num_sieve; i++) + nmod_poly_clear(sieve[i]); + nmod_poly_clear(h); + } + + return found; +} + +/* Only intended for GF(2), so we only bother with coefficients 1 */ +int +nmod_poly_irreducible_pentanomial(nmod_poly_t res, ulong n) +{ + ulong p = res->mod.n; + nmod_t mod = res->mod; + ulong k, l, m, a; + slong i; + int want_pruning; + int found = 0; + nmod_poly_t sieve[10], h; + slong num_sieve = 0; + ulong Bcoeffs[5]; + slong Bexps[5]; + ulong sieve_limit; + + if (n <= 3) + return 0; + + sieve_limit = trinomial_sieve_limit(p, n); + want_pruning = (sieve_limit > 0); + + if (want_pruning) + { + nmod_poly_init_mod(h, mod); + + for (a = 2; a <= sieve_limit && a < n; a++) + { + nmod_poly_init_mod(sieve[num_sieve], mod); + nmod_poly_product_all_irreducibles_deg(sieve[num_sieve], a); + num_sieve++; + } + } + + Bcoeffs[0] = 1; + Bcoeffs[1] = 1; + Bcoeffs[2] = 1; + Bcoeffs[3] = 1; + Bcoeffs[4] = 1; + + /* Iterate over all x^n + x^k + x^l + x^m + 1 */ + for (k = 3; k < n; k++) + { + for (l = 2; l < k; l++) + { + for (m = 1; m < l; m++) + { + /* Perfect square */ + if (p == 2 && n % 2 == 0 && k % 2 == 0 && l % 2 == 0 && m % 2 == 0) + goto next_pentanomial; + + nmod_poly_zero(res); + nmod_poly_set_coeff_ui(res, n, 1); + nmod_poly_set_coeff_ui(res, k, 1); + nmod_poly_set_coeff_ui(res, l, 1); + nmod_poly_set_coeff_ui(res, m, 1); + nmod_poly_set_coeff_ui(res, 0, 1); + + for (i = 0; i < num_sieve; i++) + { + Bexps[0] = 0; + Bexps[1] = m; + Bexps[2] = l; + Bexps[3] = k; + Bexps[4] = n; + + nmod_poly_gcd_with_sparse(h, sieve[i], res, Bcoeffs, Bexps, 5); + if (!nmod_poly_is_one(h)) + goto next_pentanomial; + } + + if (nmod_poly_is_irreducible_ddf(res)) + { + found = 1; + goto cleanup; + } + + next_pentanomial: + continue; + } + } + } + +cleanup: + if (want_pruning) + { + for (i = 0; i < num_sieve; i++) + nmod_poly_clear(sieve[i]); + nmod_poly_clear(h); + } + + return found; +} + +void +nmod_poly_minimal_irreducible(nmod_poly_t res, ulong n) +{ + FLINT_ASSERT(n != 0); + + if (n == 1) + { + nmod_poly_zero(res); + nmod_poly_set_coeff_ui(res, 1, 1); + return; + } + + if (nmod_poly_irreducible_binomial(res, n)) + return; + + if (nmod_poly_irreducible_trinomial(res, n)) + return; + + if (nmod_poly_irreducible_tetranomial(res, n)) + return; + + if (nmod_poly_irreducible_pentanomial(res, n)) + return; + + /* The conjecture is that a tetranomial (p != 2) or pentanomial (p = 2) + always works. */ + flint_throw(FLINT_ERROR, "Failed to construct minimal irreducible polynomial (n = %wu, p = %wu)", n, res->mod.n); +} + diff --git a/src/nmod_poly/test/main.c b/src/nmod_poly/test/main.c index 7146e0b8db..c6e2ab87d1 100644 --- a/src/nmod_poly/test/main.c +++ b/src/nmod_poly/test/main.c @@ -75,6 +75,7 @@ #include "t-invsqrt_series.c" #include "t-log_series.c" #include "t-make_monic.c" +#include "t-minimal_irreducible.c" #include "t-mul.c" #include "t-mul_classical.c" #include "t-mulhigh.c" @@ -199,6 +200,7 @@ test_struct tests[] = TEST_FUNCTION(nmod_poly_invsqrt_series), TEST_FUNCTION(nmod_poly_log_series), TEST_FUNCTION(nmod_poly_make_monic), + TEST_FUNCTION(nmod_poly_minimal_irreducible), TEST_FUNCTION(nmod_poly_mul), TEST_FUNCTION(nmod_poly_mul_classical), TEST_FUNCTION(nmod_poly_mulhigh), diff --git a/src/nmod_poly/test/t-minimal_irreducible.c b/src/nmod_poly/test/t-minimal_irreducible.c new file mode 100644 index 0000000000..0d5f06c5d5 --- /dev/null +++ b/src/nmod_poly/test/t-minimal_irreducible.c @@ -0,0 +1,748 @@ +/* + Copyright (C) 2025 Fredrik Johansson + + This file is part of FLINT. + + FLINT is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. See . +*/ + +#include "test_helpers.h" +#include "ulong_extras.h" +#include "nmod_poly.h" +#include "nmod_poly_factor.h" + +int +nmod_poly_irreducible_binomial_naive(nmod_poly_t res, ulong n) +{ + ulong p = res->mod.n; + ulong a; + + if (n == 0) + return 0; + + nmod_poly_zero(res); + nmod_poly_set_coeff_ui(res, n, 1); + + for (a = 1; a < p; a++) + { + if (a > 5) /* fine for the test code */ + return 0; + + nmod_poly_set_coeff_ui(res, 0, a); + if (nmod_poly_is_irreducible(res)) + return 1; + } + + return 0; +} + +int +nmod_poly_irreducible_trinomial_naive(nmod_poly_t res, ulong n) +{ + ulong p = res->mod.n; + ulong a, b, k; + + if (n <= 1) + return 0; + + /* Iterate over all x^n + a x^k + b */ + for (a = 1; a < p; a++) + { + for (b = 1; b < p; b++) + { + for (k = 1; k < n; k++) + { + if (k > 2 * n) + continue; + + nmod_poly_zero(res); + nmod_poly_set_coeff_ui(res, n, 1); + nmod_poly_set_coeff_ui(res, 0, b); + nmod_poly_set_coeff_ui(res, k, a); + if (nmod_poly_is_irreducible(res)) + return 1; + } + } + } + + return 0; +} + +int +nmod_poly_irreducible_tetranomial_naive(nmod_poly_t res, ulong n) +{ + ulong p = res->mod.n; + ulong a, b, c, k, l; + + if (n <= 2) + return 0; + + /* Iterate over all x^n + a x^k + b x^l + c */ + for (a = 1; a < p; a++) + { + for (b = 1; b < p; b++) + { + for (c = 1; c < p; c++) + { + for (k = 2; k < n; k++) + { + for (l = 1; l < k; l++) + { + nmod_poly_zero(res); + nmod_poly_set_coeff_ui(res, n, 1); + nmod_poly_set_coeff_ui(res, k, a); + nmod_poly_set_coeff_ui(res, l, b); + nmod_poly_set_coeff_ui(res, 0, c); + if (nmod_poly_is_irreducible(res)) + return 1; + } + } + } + } + } + + return 0; +} + +int +nmod_poly_irreducible_pentanomial_naive(nmod_poly_t res, ulong n) +{ + ulong k, l, m; + + if (n <= 3) + return 0; + + /* Iterate over all x^n + x^k + x^l + x^m + 1 */ + for (k = 3; k < n; k++) + { + for (l = 2; l < k; l++) + { + for (m = 1; m < l; m++) + { + nmod_poly_zero(res); + nmod_poly_set_coeff_ui(res, n, 1); + nmod_poly_set_coeff_ui(res, k, 1); + nmod_poly_set_coeff_ui(res, l, 1); + nmod_poly_set_coeff_ui(res, m, 1); + nmod_poly_set_coeff_ui(res, 0, 1); + if (nmod_poly_is_irreducible(res)) + return 1; + } + } + } + + return 0; +} + +int nmod_poly_irreducible_binomial(nmod_poly_t res, ulong n); +int nmod_poly_irreducible_trinomial(nmod_poly_t res, ulong n); +int nmod_poly_irreducible_tetranomial(nmod_poly_t res, ulong n); +int nmod_poly_irreducible_pentanomial(nmod_poly_t res, ulong n); + +int nontrinomials_tab[][3] = { + {2, 1, 1}, + {2, 8, 5}, + {2, 13, 5}, + {2, 16, 5}, + {2, 19, 5}, + {2, 24, 5}, + {2, 26, 5}, + {2, 27, 5}, + {2, 32, 5}, + {2, 37, 5}, + {2, 38, 5}, + {2, 40, 5}, + {2, 43, 5}, + {2, 45, 5}, + {2, 48, 5}, + {2, 50, 5}, + {2, 51, 5}, + {2, 53, 5}, + {2, 56, 5}, + {2, 59, 5}, + {2, 61, 5}, + {2, 64, 5}, + {2, 67, 5}, + {2, 69, 5}, + {2, 70, 5}, + {2, 72, 5}, + {2, 75, 5}, + {2, 77, 5}, + {2, 78, 5}, + {2, 80, 5}, + {2, 82, 5}, + {2, 83, 5}, + {2, 85, 5}, + {2, 88, 5}, + {2, 91, 5}, + {2, 96, 5}, + {2, 99, 5}, + {2, 101, 5}, + {2, 104, 5}, + {2, 107, 5}, + {2, 109, 5}, + {2, 112, 5}, + {2, 114, 5}, + {2, 115, 5}, + {2, 116, 5}, + {2, 117, 5}, + {2, 120, 5}, + {2, 122, 5}, + {2, 125, 5}, + {2, 128, 5}, + {2, 131, 5}, + {2, 133, 5}, + {2, 136, 5}, + {2, 138, 5}, + {2, 139, 5}, + {2, 141, 5}, + {2, 143, 5}, + {2, 144, 5}, + {2, 149, 5}, + {2, 152, 5}, + {2, 157, 5}, + {2, 158, 5}, + {2, 160, 5}, + {2, 163, 5}, + {2, 164, 5}, + {2, 165, 5}, + {2, 168, 5}, + {2, 171, 5}, + {2, 173, 5}, + {2, 176, 5}, + {2, 179, 5}, + {2, 181, 5}, + {2, 184, 5}, + {2, 187, 5}, + {2, 188, 5}, + {2, 189, 5}, + {2, 190, 5}, + {2, 192, 5}, + {2, 195, 5}, + {2, 197, 5}, + {2, 200, 5}, + {2, 203, 5}, + {2, 205, 5}, + {2, 206, 5}, + {2, 208, 5}, + {2, 211, 5}, + {2, 213, 5}, + {2, 216, 5}, + {2, 219, 5}, + {2, 221, 5}, + {2, 222, 5}, + {2, 224, 5}, + {2, 226, 5}, + {2, 227, 5}, + {2, 229, 5}, + {2, 230, 5}, + {2, 232, 5}, + {2, 235, 5}, + {2, 237, 5}, + {2, 240, 5}, + {2, 243, 5}, + {2, 245, 5}, + {2, 246, 5}, + {2, 248, 5}, + {2, 251, 5}, + {2, 254, 5}, + {2, 256, 5}, + {2, 259, 5}, + {2, 261, 5}, + {2, 262, 5}, + {2, 264, 5}, + {2, 267, 5}, + {2, 269, 5}, + {2, 272, 5}, + {2, 275, 5}, + {2, 277, 5}, + {2, 280, 5}, + {2, 283, 5}, + {2, 285, 5}, + {2, 288, 5}, + {2, 290, 5}, + {2, 291, 5}, + {2, 293, 5}, + {2, 296, 5}, + {2, 298, 5}, + {2, 299, 5}, + {2, 301, 5}, + {2, 304, 5}, + {2, 306, 5}, + {2, 307, 5}, + {2, 309, 5}, + {2, 311, 5}, + {2, 312, 5}, + {2, 315, 5}, + {2, 317, 5}, + {2, 320, 5}, + {2, 323, 5}, + {2, 325, 5}, + {2, 326, 5}, + {2, 328, 5}, + {2, 331, 5}, + {2, 334, 5}, + {2, 335, 5}, + {2, 336, 5}, + {2, 338, 5}, + {2, 339, 5}, + {2, 341, 5}, + {2, 344, 5}, + {2, 347, 5}, + {2, 349, 5}, + {2, 352, 5}, + {2, 355, 5}, + {2, 356, 5}, + {2, 357, 5}, + {2, 360, 5}, + {2, 361, 5}, + {2, 363, 5}, + {2, 365, 5}, + {2, 368, 5}, + {2, 371, 5}, + {2, 373, 5}, + {2, 374, 5}, + {2, 376, 5}, + {2, 379, 5}, + {2, 381, 5}, + {2, 384, 5}, + {2, 387, 5}, + {2, 389, 5}, + {2, 392, 5}, + {2, 395, 5}, + {2, 397, 5}, + {2, 398, 5}, + {2, 400, 5}, + {2, 403, 5}, + {2, 405, 5}, + {2, 408, 5}, + {2, 410, 5}, + {2, 411, 5}, + {2, 413, 5}, + {2, 416, 5}, + {2, 419, 5}, + {2, 421, 5}, + {2, 424, 5}, + {2, 427, 5}, + {2, 429, 5}, + {2, 430, 5}, + {2, 432, 5}, + {2, 434, 5}, + {2, 435, 5}, + {2, 437, 5}, + {2, 440, 5}, + {2, 442, 5}, + {2, 443, 5}, + {2, 445, 5}, + {2, 448, 5}, + {2, 451, 5}, + {2, 452, 5}, + {2, 453, 5}, + {2, 454, 5}, + {2, 456, 5}, + {2, 459, 5}, + {2, 461, 5}, + {2, 464, 5}, + {2, 466, 5}, + {2, 467, 5}, + {2, 469, 5}, + {2, 472, 5}, + {2, 475, 5}, + {2, 477, 5}, + {2, 480, 5}, + {2, 482, 5}, + {2, 483, 5}, + {2, 485, 5}, + {2, 488, 5}, + {2, 491, 5}, + {2, 493, 5}, + {2, 496, 5}, + {2, 499, 5}, + {3, 1, 1}, + {3, 2, 2}, + {3, 49, 4}, + {3, 57, 4}, + {3, 65, 4}, + {3, 68, 4}, + {3, 75, 4}, + {3, 98, 4}, + {3, 105, 4}, + {3, 123, 4}, + {3, 129, 4}, + {3, 130, 4}, + {3, 132, 4}, + {3, 149, 4}, + {3, 161, 4}, + {3, 175, 4}, + {3, 189, 4}, + {3, 197, 4}, + {3, 207, 4}, + {3, 212, 4}, + {3, 213, 4}, + {3, 221, 4}, + {3, 223, 4}, + {3, 231, 4}, + {3, 233, 4}, + {3, 264, 4}, + {3, 267, 4}, + {3, 276, 4}, + {3, 281, 4}, + {3, 292, 4}, + {3, 297, 4}, + {3, 298, 4}, + {3, 303, 4}, + {3, 309, 4}, + {3, 311, 4}, + {3, 319, 4}, + {3, 332, 4}, + {3, 343, 4}, + {3, 391, 4}, + {3, 394, 4}, + {3, 397, 4}, + {3, 401, 4}, + {3, 404, 4}, + {3, 405, 4}, + {3, 411, 4}, + {3, 415, 4}, + {3, 426, 4}, + {3, 435, 4}, + {3, 436, 4}, + {3, 437, 4}, + {3, 439, 4}, + {3, 441, 4}, + {3, 442, 4}, + {3, 453, 4}, + {3, 459, 4}, + {3, 463, 4}, + {3, 489, 4}, + {3, 497, 4}, + {5, 1, 1}, + {5, 2, 2}, + {5, 4, 2}, + {5, 8, 2}, + {5, 16, 2}, + {5, 32, 2}, + {5, 35, 4}, + {5, 64, 2}, + {5, 70, 4}, + {5, 123, 4}, + {5, 125, 4}, + {5, 128, 2}, + {5, 140, 4}, + {5, 181, 4}, + {5, 191, 4}, + {5, 209, 4}, + {5, 213, 4}, + {5, 219, 4}, + {5, 237, 4}, + {5, 249, 4}, + {5, 250, 4}, + {5, 253, 4}, + {5, 256, 2}, + {5, 265, 4}, + {5, 273, 4}, + {5, 280, 4}, + {5, 285, 4}, + {5, 307, 4}, + {5, 319, 4}, + {5, 345, 4}, + {5, 351, 4}, + {5, 362, 4}, + {5, 375, 4}, + {5, 382, 4}, + {5, 407, 4}, + {5, 413, 4}, + {5, 415, 4}, + {5, 418, 4}, + {5, 421, 4}, + {5, 429, 4}, + {5, 433, 4}, + {5, 441, 4}, + {5, 445, 4}, + {5, 457, 4}, + {5, 461, 4}, + {5, 498, 4}, + {5, 500, 4}, + {7, 1, 1}, + {7, 2, 2}, + {7, 3, 2}, + {7, 6, 2}, + {7, 9, 2}, + {7, 18, 2}, + {7, 27, 2}, + {7, 54, 2}, + {7, 81, 2}, + {7, 124, 4}, + {7, 162, 2}, + {7, 163, 4}, + {7, 243, 2}, + {7, 268, 4}, + {7, 301, 4}, + {7, 359, 4}, + {7, 364, 4}, + {7, 372, 4}, + {7, 385, 4}, + {7, 389, 4}, + {7, 407, 4}, + {7, 431, 4}, + {7, 476, 4}, + {7, 486, 2}, + {7, 489, 4}, + {11, 1, 1}, + {11, 2, 2}, + {11, 5, 2}, + {11, 10, 2}, + {11, 25, 2}, + {11, 50, 2}, + {11, 125, 2}, + {11, 219, 4}, + {11, 250, 2}, + {11, 291, 4}, + {11, 467, 4}, + {13, 1, 1}, + {13, 2, 2}, + {13, 3, 2}, + {13, 4, 2}, + {13, 6, 2}, + {13, 8, 2}, + {13, 9, 2}, + {13, 12, 2}, + {13, 16, 2}, + {13, 18, 2}, + {13, 24, 2}, + {13, 27, 2}, + {13, 32, 2}, + {13, 36, 2}, + {13, 48, 2}, + {13, 54, 2}, + {13, 64, 2}, + {13, 72, 2}, + {13, 81, 2}, + {13, 96, 2}, + {13, 108, 2}, + {13, 128, 2}, + {13, 144, 2}, + {13, 162, 2}, + {13, 192, 2}, + {13, 216, 2}, + {13, 243, 2}, + {13, 256, 2}, + {13, 288, 2}, + {13, 324, 2}, + {13, 384, 2}, + {13, 432, 2}, + {13, 486, 2}, + {17, 1, 1}, + {17, 2, 2}, + {17, 4, 2}, + {17, 8, 2}, + {17, 16, 2}, + {17, 32, 2}, + {17, 64, 2}, + {17, 128, 2}, + {17, 231, 4}, + {17, 256, 2}, + {17, 375, 4}, + {19, 1, 1}, + {19, 2, 2}, + {19, 3, 2}, + {19, 6, 2}, + {19, 9, 2}, + {19, 18, 2}, + {19, 27, 2}, + {19, 54, 2}, + {19, 81, 2}, + {19, 162, 2}, + {19, 243, 2}, + {257, 1, 1}, + {257, 2, 2}, + {257, 4, 2}, + {257, 8, 2}, + {257, 16, 2}, + {257, 32, 2}, + {257, 64, 2}, + {367, 1, 1}, + {367, 2, 2}, + {367, 3, 2}, + {367, 6, 2}, + {367, 9, 2}, + {367, 18, 2}, + {367, 27, 2}, + {367, 54, 2}, + {367, 61, 2}, + {367, 81, 2}, + {1709, 1, 1}, + {1709, 2, 2}, + {1709, 4, 2}, + {1709, 7, 2}, + {1709, 8, 2}, + {1709, 14, 2}, + {1709, 16, 2}, + {1709, 28, 2}, + {1709, 32, 2}, + {1709, 49, 2}, + {1709, 56, 2}, + {1709, 61, 2}, + {1709, 64, 2}, + {1709, 98, 2}, + {13037, 1, 1}, + {13037, 2, 2}, + {13037, 4, 2}, + {13037, 8, 2}, + {13037, 16, 2}, + {13037, 32, 2}, + {13037, 64, 2}, + {15161, 1, 1}, + {15161, 2, 2}, + {15161, 4, 2}, + {15161, 5, 2}, + {15161, 8, 2}, + {15161, 10, 2}, + {15161, 16, 2}, + {15161, 20, 2}, + {15161, 25, 2}, + {15161, 32, 2}, + {15161, 40, 2}, + {15161, 50, 2}, + {15161, 64, 2}, + {15161, 80, 2}, + {15161, 100, 2}, + {69337, 1, 1}, + {69337, 2, 2}, + {69337, 3, 2}, + {69337, 4, 2}, + {69337, 6, 2}, + {69337, 8, 2}, + {69337, 9, 2}, + {69337, 12, 2}, + {69337, 16, 2}, + {69337, 18, 2}, + {69337, 24, 2}, + {69337, 27, 2}, + {69337, 32, 2}, + {69337, 36, 2}, + {69337, 48, 2}, + {69337, 54, 2}, + {69337, 64, 2}, + {69337, 72, 2}, + {69337, 81, 2}, + {69337, 96, 2}, + {438467, 1, 1}, + {438467, 2, 2}, + {438467, 7, 2}, + {438467, 14, 2}, + {438467, 49, 2}, + {438467, 98, 2}, + {0, 0, 0}, +}; + + +TEST_FUNCTION_START(nmod_poly_minimal_irreducible, state) +{ + ulong p, n, nmax; + slong c1, c2, j, i, pi; + nmod_poly_t f, g; + int res1, res2; + + ulong ps[] = { 2, 3, 5, 7, 11, 17, 19, 29, 31, 37, 10007, UWORD(2147483659), +#if FLINT_BITS == 64 + UWORD(9223372036854775837), +#endif + 0 }; + + for (i = 0; nontrinomials_tab[i][0] != 0; i++) + { + p = nontrinomials_tab[i][0]; + n = nontrinomials_tab[i][1]; + c2 = nontrinomials_tab[i][2]; + + if (n > FLINT_MAX(50, flint_test_multiplier() * 10)) + continue; + + /* flint_printf("%wu %wu\n", p, n); */ + + nmod_poly_init(f, p); + nmod_poly_minimal_irreducible(f, n); + c1 = 0; + for (j = 0; j < f->length; j++) + c1 += (f->coeffs[j] != 0); + + if (c1 != c2) + { + flint_printf("FAIL (wrong number of terms): p = %wu n = %wu c1 = %wd c2 = %wd f = %{nmod_poly}\n", + p, n, c1, c2, f); + flint_abort(); + } + + nmod_poly_clear(f); + } + + for (pi = 0; (p = ps[pi]) != 0; pi++) + { + nmod_poly_init(f, p); + nmod_poly_init(g, p); + + if (p == 2) + nmax = 200; + else if (p <= 5) + nmax = 100; + else + nmax = 10; + + nmax = FLINT_MIN(nmax, FLINT_MAX(50, flint_test_multiplier() * 10)); + + /* flint_printf("%wu\n", p); */ + + for (n = 0; n <= nmax; n++) + { + /* flint_printf("n = %wd\n", n); */ + + res1 = nmod_poly_irreducible_binomial(f, n); + res2 = nmod_poly_irreducible_binomial_naive(g, n); + + if (res1 != res2 || (res1 && !nmod_poly_equal(f, g))) + { + flint_printf("FAIL (binomial): p = %wu n = %wu res1 = %d res2 = %d\n", p, n, res1, res2); + flint_abort(); + } + + res1 = nmod_poly_irreducible_trinomial(f, n); + res2 = nmod_poly_irreducible_trinomial_naive(g, n); + + if (res1 != res2) + { + flint_printf("FAIL (trinomial): p = %wu n = %wu res1 = %d res2 = %d\n", p, n, res1, res2); + flint_abort(); + } + + res1 = nmod_poly_irreducible_tetranomial(f, n); + res2 = nmod_poly_irreducible_tetranomial_naive(g, n); + + if (res1 != res2) + { + flint_printf("FAIL (tetranomial): p = %wu n = %wu res1 = %d res2 = %d\n", p, n, res1, res2); + flint_abort(); + } + + if (p == 2) + { + res1 = nmod_poly_irreducible_pentanomial(f, n); + res2 = nmod_poly_irreducible_pentanomial_naive(g, n); + + if (res1 != res2) + { + flint_printf("FAIL (pentanomial): p = %wu n = %wu res1 = %d res2 = %d\n", p, n, res1, res2); + flint_abort(); + } + } + } + + nmod_poly_clear(f); + nmod_poly_clear(g); + } + + TEST_FUNCTION_END(state); +}