Ulli/flint2
Ulli/flint2
1
0
Fork 0
Code Issues Pull requests Packages Projects Releases Wiki Activity

add fmpz_mpoly_resultant and discriminant

Browse source
This commit is contained in:
Daniel Schultz 2021-06-01 22:21:25 +02:00
parent bc5616f099
commit 4347e7e473
21 changed files with 1076 additions and 134 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

8
doc/source/fmpz_mpoly.rst
View file

@ -611,6 +611,14 @@ Greatest Common Divisor
Try to set ``G`` to the GCD of ``A`` and ``B`` using various algorithms.
.. function:: int fmpz_mpoly_resultant(fmpz_mpoly_t R, const fmpz_mpoly_t A, const fmpz_mpoly_t B, slong var, const fmpz_mpoly_ctx_t ctx)
Try to set ``R`` to the resultant of ``A`` and ``B`` with respect to the variable of index ``var``.
.. function:: int fmpz_mpoly_discriminant(fmpz_mpoly_t D, const fmpz_mpoly_t A, slong var, const fmpz_mpoly_ctx_t ctx)
Try to set ``D`` to the discriminant of ``A`` and ``B`` with respect to the variable of index ``var``.
Square Root
--------------------------------------------------------------------------------

7
fmpz_mod_mpoly.h
View file

@ -726,11 +726,8 @@ FLINT_DLL void fmpz_mod_mpoly_divrem_ideal(fmpz_mod_mpoly_struct ** Q,
const fmpz_mod_mpoly_ctx_t ctx);
FMPZ_MOD_MPOLY_INLINE
void fmpz_mod_mpoly_divexact(
fmpz_mod_mpoly_t Q,
const fmpz_mod_mpoly_t A,
const fmpz_mod_mpoly_t B,
const fmpz_mod_mpoly_ctx_t ctx)
void fmpz_mod_mpoly_divexact(fmpz_mod_mpoly_t Q, const fmpz_mod_mpoly_t A,
const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_ctx_t ctx)
{
if (fmpz_mod_mpoly_divides(Q, A, B, ctx))
return;

3
fmpz_mod_mpoly/get_set_is_fmpz_mod_poly.c
View file

@ -10,9 +10,6 @@
*/
#include "fmpz_mod_mpoly_factor.h"
#include "nmod_mpoly_factor.h"
int fmpz_mod_mpoly_is_fmpz_mod_poly(
const fmpz_mod_mpoly_t A,

6
fmpz_mod_mpoly/univar.c
View file

@ -118,8 +118,6 @@ void fmpz_mod_mpoly_univar_set_coeff_ui(
return;
}
void fmpz_mod_mpoly_univar_assert_canonical(
fmpz_mod_mpoly_univar_t A,
const fmpz_mod_mpoly_ctx_t ctx)
@ -569,7 +567,7 @@ void fmpz_mod_mpoly_univar_set(
#define COEFF(A, i) ((void*)(A->coeffs + (i)*R->elem_size))
void mpoly_univar_set_fmpz_mod_mpoly_univar(
static void mpoly_univar_set_fmpz_mod_mpoly_univar(
mpoly_univar_t A,
mpoly_void_ring_t R,
const fmpz_mod_mpoly_univar_t B,
@ -587,7 +585,7 @@ void mpoly_univar_set_fmpz_mod_mpoly_univar(
}
}
void mpoly_univar_swap_fmpz_mod_mpoly_univar(
static void mpoly_univar_swap_fmpz_mod_mpoly_univar(
mpoly_univar_t A,
mpoly_void_ring_t R,
fmpz_mod_mpoly_univar_t B,

51
fmpz_mpoly.h
View file

@ -467,6 +467,20 @@ FLINT_DLL void fmpz_mpoly_get_coeff_vars_ui(fmpz_mpoly_t C,
const fmpz_mpoly_t A, const slong * vars, const ulong * exps,
slong length, const fmpz_mpoly_ctx_t ctx);
/* conversion ****************************************************************/
FLINT_DLL int fmpz_mpoly_is_fmpz_poly(const fmpz_mpoly_t A, slong var,
const fmpz_mpoly_ctx_t ctx);
FLINT_DLL int fmpz_mpoly_get_fmpz_poly(fmpz_poly_t A, const fmpz_mpoly_t B,
slong var, const fmpz_mpoly_ctx_t ctx);
FLINT_DLL void _fmpz_mpoly_set_fmpz_poly(fmpz_mpoly_t A, flint_bitcnt_t Abits,
const fmpz * Bcoeffs, slong Blen, slong var, const fmpz_mpoly_ctx_t ctx);
FLINT_DLL void fmpz_mpoly_set_fmpz_poly(fmpz_mpoly_t A, const fmpz_poly_t B,
slong v, const fmpz_mpoly_ctx_t ctx);
/* comparison ****************************************************************/
FLINT_DLL int fmpz_mpoly_cmp(const fmpz_mpoly_t A, const fmpz_mpoly_t B,
@ -859,6 +873,15 @@ FLINT_DLL void fmpz_mpoly_quasidivrem_ideal(fmpz_t scale,
fmpz_mpoly_struct ** Q, fmpz_mpoly_t R, const fmpz_mpoly_t A,
fmpz_mpoly_struct * const * B, slong len, const fmpz_mpoly_ctx_t ctx);
FMPZ_MPOLY_INLINE
void fmpz_mpoly_divexact(fmpz_mpoly_t Q, const fmpz_mpoly_t A,
const fmpz_mpoly_t B, const fmpz_mpoly_ctx_t ctx)
{
if (fmpz_mpoly_divides(Q, A, B, ctx))
return;
flint_throw(FLINT_ERROR, "fmpz_mpoly_divexact: nonexact division");
}
FLINT_DLL slong _fmpz_mpoly_div_monagan_pearce(fmpz ** polyq,
ulong ** expq, slong * allocq, const fmpz * poly2,
@ -997,10 +1020,19 @@ FLINT_DLL void fmpz_mpoly_univar_print_pretty(const fmpz_mpoly_univar_t A,
FLINT_DLL void fmpz_mpoly_univar_assert_canonical(fmpz_mpoly_univar_t A,
const fmpz_mpoly_ctx_t ctx);
FMPZ_MPOLY_INLINE
void fmpz_mpoly_univar_zero(fmpz_mpoly_univar_t A, const fmpz_mpoly_ctx_t ctx)
{
A->length = 0;
}
FLINT_DLL void fmpz_mpoly_univar_set_coeff_ui(fmpz_mpoly_univar_t A,
ulong e, const fmpz_mpoly_t c, const fmpz_mpoly_ctx_t ctx);
FLINT_DLL void fmpz_mpoly_to_univar(fmpz_mpoly_univar_t A,
const fmpz_mpoly_t B, slong var, const fmpz_mpoly_ctx_t ctx);
FLINT_DLL void fmpz_mpoly_from_univar_bits(fmpz_mpoly_t A, flint_bitcnt_t Abits,
FLINT_DLL void _fmpz_mpoly_from_univar(fmpz_mpoly_t A, flint_bitcnt_t Abits,
const fmpz_mpoly_univar_t B, slong var, const fmpz_mpoly_ctx_t ctx);
FLINT_DLL void fmpz_mpoly_from_univar(fmpz_mpoly_t A,
@ -1053,6 +1085,23 @@ void fmpz_mpoly_univar_swap_term_coeff(fmpz_mpoly_t c,
fmpz_mpoly_swap(c, A->coeffs + i, ctx);
}
FLINT_DLL int fmpz_mpoly_univar_pseudo_gcd(fmpz_mpoly_univar_t gx,
const fmpz_mpoly_univar_t ax, const fmpz_mpoly_univar_t bx,
const fmpz_mpoly_ctx_t ctx);
FLINT_DLL int fmpz_mpoly_univar_resultant(fmpz_mpoly_t d,
const fmpz_mpoly_univar_t ax, const fmpz_mpoly_univar_t bx,
const fmpz_mpoly_ctx_t ctx);
FLINT_DLL int fmpz_mpoly_univar_discriminant(fmpz_mpoly_t d,
const fmpz_mpoly_univar_t fx, const fmpz_mpoly_ctx_t ctx);
FLINT_DLL int fmpz_mpoly_resultant(fmpz_mpoly_t R, const fmpz_mpoly_t A,
const fmpz_mpoly_t B, slong var, const fmpz_mpoly_ctx_t ctx);
FLINT_DLL int fmpz_mpoly_discriminant(fmpz_mpoly_t R, const fmpz_mpoly_t A,
slong var, const fmpz_mpoly_ctx_t ctx);
/******************************************************************************
Internal functions (guaranteed to change without notice)

30
fmpz_mpoly/discriminant.c Normal file
View file

@ -0,0 +1,30 @@
/*
Copyright (C) 2021 Daniel Schultz
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 2.1 of the License, or
(at your option) any later version. See <https://www.gnu.org/licenses/>.
*/
#include "fmpz_mpoly.h"
int fmpz_mpoly_discriminant(fmpz_mpoly_t R, const fmpz_mpoly_t A,
slong var, const fmpz_mpoly_ctx_t ctx)
{
int success;
fmpz_mpoly_univar_t Ax;
fmpz_mpoly_univar_init(Ax, ctx);
fmpz_mpoly_to_univar(Ax, A, var, ctx);
success = fmpz_mpoly_univar_discriminant(R, Ax, ctx);
fmpz_mpoly_univar_clear(Ax, ctx);
return success;
}

43
fmpz_mpoly_factor/get_set_is_fmpz_poly.c → fmpz_mpoly/get_set_is_fmpz_poly.c
View file

@ -9,9 +9,7 @@
(at your option) any later version. See <https://www.gnu.org/licenses/>.
*/
#include "fmpz_mpoly_factor.h"
#include "nmod_mpoly_factor.h"
#include "fmpz_mpoly.h"
int fmpz_mpoly_is_fmpz_poly(
const fmpz_mpoly_t A,
@ -21,7 +19,6 @@ int fmpz_mpoly_is_fmpz_poly(
return mpoly_is_poly(A->exps, A->length, A->bits, var, ctx->minfo);
}
int fmpz_mpoly_get_fmpz_poly(
fmpz_poly_t A,
const fmpz_mpoly_t B,
@ -77,42 +74,6 @@ int fmpz_mpoly_get_fmpz_poly(
}
}
void fmpz_mpoly_fit_length_set_bits(
fmpz_mpoly_t A,
slong len,
flint_bitcnt_t bits,
const fmpz_mpoly_ctx_t ctx)
{
slong N = mpoly_words_per_exp(bits, ctx->minfo);
slong new_alloc;
FLINT_ASSERT(len >= 0);
if (A->alloc < len)
{
new_alloc = FLINT_MAX(len, 2*A->alloc);
if (A->alloc > 0)
{
A->coeffs = (fmpz *) flint_realloc(A->coeffs, new_alloc*sizeof(fmpz));
A->exps = (ulong *) flint_realloc(A->exps, new_alloc*N*sizeof(ulong));
memset(A->coeffs + A->alloc, 0, (new_alloc - A->alloc)*sizeof(fmpz));
}
else
{
A->coeffs = (fmpz *) flint_calloc(new_alloc, sizeof(fmpz));
A->exps = (ulong *) flint_malloc(new_alloc*N*sizeof(ulong));
}
A->alloc = new_alloc;
}
else if (A->bits < bits)
{
if (A->alloc > 0)
A->exps = (ulong *) flint_realloc(A->exps, A->alloc*N*sizeof(ulong));
}
A->bits = bits;
}
void _fmpz_mpoly_set_fmpz_poly(
fmpz_mpoly_t A,
flint_bitcnt_t Abits,
@ -138,7 +99,7 @@ void _fmpz_mpoly_set_fmpz_poly(
for (i = 0; i < Blen; i++)
Alen += (Bcoeffs[i] != 0);
fmpz_mpoly_fit_length_set_bits(A, Alen, Abits, ctx);
fmpz_mpoly_fit_length_reset_bits(A, Alen, Abits, ctx);
Alen = 0;
for (i = Blen - 1; i >= 0; i--)

33
fmpz_mpoly/resultant.c Normal file
View file

@ -0,0 +1,33 @@
/*
Copyright (C) 2021 Daniel Schultz
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 2.1 of the License, or
(at your option) any later version. See <https://www.gnu.org/licenses/>.
*/
#include "fmpz_mpoly.h"
int fmpz_mpoly_resultant(fmpz_mpoly_t R, const fmpz_mpoly_t A,
const fmpz_mpoly_t B, slong var, const fmpz_mpoly_ctx_t ctx)
{
int success;
fmpz_mpoly_univar_t Ax, Bx;
fmpz_mpoly_univar_init(Ax, ctx);
fmpz_mpoly_univar_init(Bx, ctx);
fmpz_mpoly_to_univar(Ax, A, var, ctx);
fmpz_mpoly_to_univar(Bx, B, var, ctx);
success = fmpz_mpoly_univar_resultant(R, Ax, Bx, ctx);
fmpz_mpoly_univar_clear(Ax, ctx);
fmpz_mpoly_univar_clear(Bx, ctx);
return success;
}

256
fmpz_mpoly/test/t-resultant_discriminant.c Normal file
View file

@ -0,0 +1,256 @@
/*
Copyright (C) 2021 Daniel Schultz
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 2.1 of the License, or
(at your option) any later version. See <https://www.gnu.org/licenses/>.
*/
#include <stdio.h>
#include <stdlib.h>
#include "fmpz_mpoly.h"
int
main(void)
{
slong i, j;
FLINT_TEST_INIT(state);
flint_printf("resultant_discriminant....");
fflush(stdout);
/* Check quadratic polynomial */
{
fmpz_mpoly_ctx_t ctx;
fmpz_mpoly_t f, d, d1;
const char * vars[] = {"x","a","b","c"};
fmpz_mpoly_ctx_init(ctx, 4, ORD_DEGLEX);
fmpz_mpoly_init(f, ctx);
fmpz_mpoly_init(d, ctx);
fmpz_mpoly_init(d1, ctx);
fmpz_mpoly_set_str_pretty(f, "a^10*x^2 + b^100*x + c^100000000000000000000", vars, ctx);
fmpz_mpoly_set_str_pretty(d1, "b^200 - 4*a^10*c^100000000000000000000", vars, ctx);
if (!fmpz_mpoly_discriminant(d, f, 0, ctx))
{
flint_printf("FAIL: could not compute quadratic discriminant\n");
flint_abort();
}
if (!fmpz_mpoly_equal(d, d1, ctx))
{
flint_printf("FAIL: Check quadratic polynomial\n");
flint_abort();
}
fmpz_mpoly_clear(f, ctx);
fmpz_mpoly_clear(d, ctx);
fmpz_mpoly_clear(d1, ctx);
fmpz_mpoly_ctx_clear(ctx);
}
/* Check univariate resultant */
for (i = 0; i < 50 * flint_test_multiplier(); i++)
{
fmpz_mpoly_ctx_t ctx;
fmpz_mpoly_t a, b, r;
fmpz_poly_t au, bu;
fmpz_t ru;
ordering_t ord;
slong nvars, len1, len2, exp_bound1, exp_bound2;
slong coeff_bits;
ord = mpoly_ordering_randtest(state);
nvars = 1;
fmpz_mpoly_ctx_init(ctx, nvars, ord);
fmpz_mpoly_init(a, ctx);
fmpz_mpoly_init(b, ctx);
fmpz_mpoly_init(r, ctx);
fmpz_poly_init(au);
fmpz_poly_init(bu);
fmpz_init(ru);
len1 = n_randint(state, 100);
len2 = n_randint(state, 100);
exp_bound1 = n_randint(state, 50) + 1;
exp_bound2 = n_randint(state, 50) + 1;
coeff_bits = n_randint(state, 4);
fmpz_mpoly_randtest_bound(a, state, len1, coeff_bits, exp_bound1, ctx);
fmpz_mpoly_randtest_bound(b, state, len2, coeff_bits, exp_bound2, ctx);
fmpz_mpoly_get_fmpz_poly(au, a, 0, ctx);
fmpz_mpoly_get_fmpz_poly(bu, b, 0, ctx);
fmpz_poly_resultant(ru, au, bu);
if (!fmpz_mpoly_resultant(r, a, b, 0, ctx) ||
!fmpz_mpoly_equal_fmpz(r, ru, ctx))
{
flint_printf("FAIL: Check univariate resultant \n");
flint_printf("i: %wd\n",i);
flint_abort();
}
fmpz_mpoly_clear(a, ctx);
fmpz_mpoly_clear(b, ctx);
fmpz_mpoly_clear(r, ctx);
fmpz_poly_clear(au);
fmpz_poly_clear(bu);
fmpz_clear(ru);
}
/* Check res(a*b,c) = res(a,c)*res(b,c) */
for (i = 0; i < 30 * flint_test_multiplier(); i++)
{
fmpz_mpoly_ctx_t ctx;
fmpz_mpoly_t a, b, c, ab, ra, rb, rab, p;
ordering_t ord;
slong nvars, len1, len2, len3, exp_bound1, exp_bound2, exp_bound3;
slong coeff_bits;
ord = mpoly_ordering_randtest(state);
nvars = n_randint(state, 3) + 1;
fmpz_mpoly_ctx_init(ctx, nvars, ord);
fmpz_mpoly_init(a, ctx);
fmpz_mpoly_init(b, ctx);
fmpz_mpoly_init(c, ctx);
fmpz_mpoly_init(ab, ctx);
fmpz_mpoly_init(ra, ctx);
fmpz_mpoly_init(rb, ctx);
fmpz_mpoly_init(rab, ctx);
fmpz_mpoly_init(p, ctx);
len1 = n_randint(state, 15);
len2 = n_randint(state, 15);
len3 = n_randint(state, 15);
exp_bound1 = n_randint(state, 5) + 1;
exp_bound2 = n_randint(state, 5) + 1;
exp_bound3 = n_randint(state, 5) + 1;
coeff_bits = n_randint(state, 3);
fmpz_mpoly_randtest_bound(a, state, len1, coeff_bits, exp_bound1, ctx);
fmpz_mpoly_randtest_bound(b, state, len2, coeff_bits, exp_bound2, ctx);
fmpz_mpoly_randtest_bound(c, state, len3, coeff_bits, exp_bound3, ctx);
for (j = 0; j < nvars; j++)
{
fmpz_mpoly_mul(ab, a, b, ctx);
if (!fmpz_mpoly_resultant(ra, a, c, j, ctx))
continue;
fmpz_mpoly_assert_canonical(ra, ctx);
if (!fmpz_mpoly_resultant(rb, b, c, j, ctx))
continue;
fmpz_mpoly_assert_canonical(rb, ctx);
if (!fmpz_mpoly_resultant(rab, ab, c, j, ctx))
continue;
fmpz_mpoly_assert_canonical(rab, ctx);
fmpz_mpoly_mul(p, ra, rb, ctx);
if (!fmpz_mpoly_equal(p,rab,ctx))
{
flint_printf("FAIL: Check res(a*b,c) = res(a,c)*res(b,c)\n");
flint_printf("i: %wd j: %wd\n",i,j);
flint_abort();
}
}
fmpz_mpoly_clear(a, ctx);
fmpz_mpoly_clear(b, ctx);
fmpz_mpoly_clear(c, ctx);
fmpz_mpoly_clear(ab, ctx);
fmpz_mpoly_clear(ra, ctx);
fmpz_mpoly_clear(rb, ctx);
fmpz_mpoly_clear(rab, ctx);
fmpz_mpoly_clear(p, ctx);
}
/* Check disc(a*b) = disc(a)*disc(b)*res(a,b)^2 */
for (i = 0; i < 30 * flint_test_multiplier(); i++)
{
fmpz_mpoly_ctx_t ctx;
fmpz_mpoly_t a, b, ab, r, da, db, dab, p;
ordering_t ord;
slong nvars, len1, len2, exp_bound1, exp_bound2;
slong coeff_bits;
ord = mpoly_ordering_randtest(state);
nvars = n_randint(state, 3) + 1;
fmpz_mpoly_ctx_init(ctx, nvars, ord);
fmpz_mpoly_init(a, ctx);
fmpz_mpoly_init(b, ctx);
fmpz_mpoly_init(ab, ctx);
fmpz_mpoly_init(da, ctx);
fmpz_mpoly_init(db, ctx);
fmpz_mpoly_init(dab, ctx);
fmpz_mpoly_init(r, ctx);
fmpz_mpoly_init(p, ctx);
len1 = n_randint(state, 15);
len2 = n_randint(state, 15);
exp_bound1 = n_randint(state, 5) + 1;
exp_bound2 = n_randint(state, 5) + 1;
coeff_bits = n_randint(state, 3);
fmpz_mpoly_randtest_bound(a, state, len1, coeff_bits, exp_bound1, ctx);
fmpz_mpoly_randtest_bound(b, state, len2, coeff_bits, exp_bound2, ctx);
for (j = 0; j < nvars; j++)
{
if (fmpz_mpoly_degree_si(a, j, ctx) < 1)
continue;
if (fmpz_mpoly_degree_si(b, j, ctx) < 1)
continue;
fmpz_mpoly_mul(ab, a, b, ctx);
if (!fmpz_mpoly_resultant(r, a, b, j, ctx))
continue;
if (!fmpz_mpoly_discriminant(da, a, j, ctx))
continue;
if (!fmpz_mpoly_discriminant(db, b, j, ctx))
continue;
if (!fmpz_mpoly_discriminant(dab, ab, j, ctx))
continue;
fmpz_mpoly_mul(p, da, db, ctx);
fmpz_mpoly_mul(p, p, r, ctx);
fmpz_mpoly_mul(p, p, r, ctx);
if (!fmpz_mpoly_equal(dab, p, ctx))
{
flint_printf("FAIL: Check disc(a*b) = disc(a)*disc(b)*res(a,b)^2\n");
flint_printf("i: %wd j: %wd\n",i,j);
flint_abort();
}
}
fmpz_mpoly_clear(a, ctx);
fmpz_mpoly_clear(b, ctx);
fmpz_mpoly_clear(ab, ctx);
fmpz_mpoly_clear(da, ctx);
fmpz_mpoly_clear(db, ctx);
fmpz_mpoly_clear(dab, ctx);
fmpz_mpoly_clear(r, ctx);
fmpz_mpoly_clear(p, ctx);
}
FLINT_TEST_CLEANUP(state);
printf("PASS\n");
return 0;
}

2
fmpz_mpoly/test/t-univar.c
View file

@ -68,7 +68,7 @@ main(void)
}
bits = mpoly_fix_bits(f->bits + n_randint(state, FLINT_BITS), ctx->minfo);
fmpz_mpoly_from_univar_bits(h, bits, fx, j, ctx);
_fmpz_mpoly_from_univar(h, bits, fx, j, ctx);
fmpz_mpoly_assert_canonical(h, ctx);
if (h->bits != bits || !fmpz_mpoly_equal(f, h, ctx))
{

180
fmpz_mpoly/test/t-univar_resultant.c Normal file
View file

@ -0,0 +1,180 @@
/*
Copyright (C) 2021 Daniel Schultz
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 2.1 of the License, or
(at your option) any later version. See <https://www.gnu.org/licenses/>.
*/
#include <stdio.h>
#include <stdlib.h>
#include "fmpz_mpoly.h"
void test_resultant(
const fmpz_mpoly_t fx,
const fmpz_mpoly_t gx,
const fmpz_mpoly_ctx_t ctx)
{
fmpz_mpoly_univar_t F, G;
fmpz_poly_t f, g;
fmpz_mpoly_t R;
fmpz_t r;
fmpz_mpoly_univar_init(F, ctx);
fmpz_mpoly_univar_init(G, ctx);
fmpz_poly_init(f);
fmpz_poly_init(g);
fmpz_mpoly_init(R, ctx);
fmpz_init(r);
fmpz_mpoly_get_fmpz_poly(f, fx, 0, ctx);
fmpz_mpoly_get_fmpz_poly(g, gx, 0, ctx);
fmpz_poly_resultant(r, f, g);
fmpz_mpoly_to_univar(F, fx, 0, ctx);
fmpz_mpoly_to_univar(G, gx, 0, ctx);
fmpz_mpoly_univar_resultant(R, F, G, ctx);
if (!fmpz_mpoly_equal_fmpz(R, r, ctx))
{
flint_printf("FAIL: Check resultant against univariate\n");
flint_printf("fx: ");
fmpz_mpoly_print_pretty(fx, NULL, ctx);
flint_printf("\n");
flint_printf("gx: ");
fmpz_mpoly_print_pretty(gx, NULL, ctx);
flint_printf("\n");
flint_printf("R: ");
fmpz_mpoly_print_pretty(R, NULL, ctx);
flint_printf("\n");
flint_printf("r: ");
fmpz_print(r);
flint_printf("\n");
flint_abort();
}
fmpz_mpoly_univar_clear(F, ctx);
fmpz_mpoly_univar_clear(G, ctx);
fmpz_poly_clear(f);
fmpz_poly_clear(g);
fmpz_mpoly_clear(R, ctx);
fmpz_clear(r);
}
int
main(void)
{
slong i, j;
FLINT_TEST_INIT(state);
flint_printf("univar_resultant....");
fflush(stdout);
{
fmpz_mpoly_ctx_t ctx;
fmpz_mpoly_t t, s;
fmpz_mpoly_univar_t f;
const char * vars[] = {"a", "b", "c", "d"};
fmpz_mpoly_ctx_init(ctx, 4, ORD_DEGLEX);
fmpz_mpoly_init(t, ctx);
fmpz_mpoly_init(s, ctx);
fmpz_mpoly_univar_init(f, ctx);
fmpz_mpoly_univar_zero(f, ctx);
fmpz_mpoly_set_str_pretty(t, "a", vars, ctx);
fmpz_mpoly_univar_set_coeff_ui(f, 1, t, ctx);
fmpz_mpoly_set_str_pretty(t, "b", vars, ctx);
fmpz_mpoly_univar_set_coeff_ui(f, 0, t, ctx);
fmpz_mpoly_set_str_pretty(t, "1", vars, ctx);
fmpz_mpoly_univar_discriminant(s, f, ctx);
if (!fmpz_mpoly_equal(s, t, ctx))
{
flint_printf("FAIL: check linear discriminant\n");
flint_abort();
}
fmpz_mpoly_univar_zero(f, ctx);
fmpz_mpoly_set_str_pretty(t, "a", vars, ctx);
fmpz_mpoly_univar_set_coeff_ui(f, 2, t, ctx);
fmpz_mpoly_set_str_pretty(t, "b", vars, ctx);
fmpz_mpoly_univar_set_coeff_ui(f, 1, t, ctx);
fmpz_mpoly_set_str_pretty(t, "c", vars, ctx);
fmpz_mpoly_univar_set_coeff_ui(f, 0, t, ctx);
fmpz_mpoly_set_str_pretty(t, "b^2-4*a*c", vars, ctx);
fmpz_mpoly_univar_discriminant(s, f, ctx);
if (!fmpz_mpoly_equal(s, t, ctx))
{
flint_printf("FAIL: check quadratic discriminant\n");
flint_abort();
}
fmpz_mpoly_univar_zero(f, ctx);
fmpz_mpoly_set_str_pretty(t, "a", vars, ctx);
fmpz_mpoly_univar_set_coeff_ui(f, 3, t, ctx);
fmpz_mpoly_set_str_pretty(t, "b", vars, ctx);
fmpz_mpoly_univar_set_coeff_ui(f, 2, t, ctx);
fmpz_mpoly_set_str_pretty(t, "c", vars, ctx);
fmpz_mpoly_univar_set_coeff_ui(f, 1, t, ctx);
fmpz_mpoly_set_str_pretty(t, "d", vars, ctx);
fmpz_mpoly_univar_set_coeff_ui(f, 0, t, ctx);
fmpz_mpoly_set_str_pretty(t, "b^2*c^2-4*a*c^3-4*b^3*d-27*a^2*d^2+18*a*b*c*d", vars, ctx);
fmpz_mpoly_univar_discriminant(s, f, ctx);
if (!fmpz_mpoly_equal(s, t, ctx))
{
flint_printf("FAIL: check cubic discriminant\n");
flint_abort();
}
fmpz_mpoly_clear(t, ctx);
fmpz_mpoly_clear(s, ctx);
fmpz_mpoly_univar_clear(f, ctx);
fmpz_mpoly_ctx_clear(ctx);
}
for (i = 0; i < 5 * flint_test_multiplier(); i++)
{
fmpz_mpoly_ctx_t ctx;
fmpz_mpoly_t f, g, t;
fmpz_mpoly_ctx_init_rand(ctx, state, 1);
fmpz_mpoly_init(f, ctx);
fmpz_mpoly_init(g, ctx);
fmpz_mpoly_init(t, ctx);
for (j = 0; j < 5; j++)
{
fmpz_mpoly_randtest_bound(f, state, 3, n_randint(state, 100) + 1, 3, ctx);
if (fmpz_mpoly_is_zero(f, ctx))
fmpz_mpoly_one(f, ctx);
fmpz_mpoly_zero(g, ctx);
while (fmpz_mpoly_degree_si(f, 0, ctx) < 100)
{
fmpz_mpoly_randtest_bound(t, state, 5, n_randint(state, 100) + 1, 10, ctx);
fmpz_mpoly_mul(t, t, f, ctx);
fmpz_mpoly_add(g, g, t, ctx);
fmpz_mpoly_swap(f, g, ctx);
}
fmpz_mpoly_randtest_bound(f, state, 20, n_randint(state, 100) + 1, 50, ctx);
fmpz_mpoly_randtest_bound(g, state, 20, n_randint(state, 100) + 1, 50, ctx);
test_resultant(f, g, ctx);
}
fmpz_mpoly_clear(f, ctx);
fmpz_mpoly_clear(g, ctx);
fmpz_mpoly_clear(t, ctx);
fmpz_mpoly_ctx_clear(ctx);
}
FLINT_TEST_CLEANUP(state);
flint_printf("PASS\n");
return 0;
}

171
fmpz_mpoly/univar.c
View file

@ -74,6 +74,58 @@ void fmpz_mpoly_univar_fit_length(fmpz_mpoly_univar_t A,
}
}
void fmpz_mpoly_univar_set_coeff_ui(
fmpz_mpoly_univar_t A,
ulong e,
const fmpz_mpoly_t c,
const fmpz_mpoly_ctx_t ctx)
{
slong i, j;
for (i = A->length; i >= 0; i--)
{
int cmp = i > 0 ? fmpz_cmp_ui(A->exps + i - 1, e) : 1;
if (cmp > 0)
{
if (fmpz_mpoly_is_zero(c, ctx))
return;
fmpz_mpoly_univar_fit_length(A, A->length + 1, ctx);
for (j = A->length; j > i; j--)
{
fmpz_mpoly_swap(A->coeffs + j, A->coeffs + j + 1, ctx);
fmpz_swap(A->exps + j, A->exps + j + 1);
}
A->length++;
fmpz_set_ui(A->exps + i, e);
fmpz_mpoly_set(A->coeffs + i, c, ctx);
return;
}
else if (cmp == 0)
{
fmpz_mpoly_set(A->coeffs + i, c, ctx);
if (!fmpz_mpoly_is_zero(A->coeffs + i, ctx))
return;
A->length--;
for (j = i; j < A->length; j++)
{
fmpz_mpoly_swap(A->coeffs + j, A->coeffs + j + 1, ctx);
fmpz_swap(A->exps + j, A->exps + j + 1);
}
}
}
FLINT_ASSERT(0 && "unreachable");
return;
}
void fmpz_mpoly_univar_assert_canonical(fmpz_mpoly_univar_t A, const fmpz_mpoly_ctx_t ctx)
{
slong i;
@ -290,7 +342,7 @@ void fmpz_mpoly_to_univar(fmpz_mpoly_univar_t A, const fmpz_mpoly_t B,
The assertion x->next == NULL would need to be replaced by a loop.
Other asserts would need to be removed as well.
*/
void fmpz_mpoly_from_univar_bits(fmpz_mpoly_t A, flint_bitcnt_t Abits,
void _fmpz_mpoly_from_univar(fmpz_mpoly_t A, flint_bitcnt_t Abits,
const fmpz_mpoly_univar_t B, slong var, const fmpz_mpoly_ctx_t ctx)
{
slong N = mpoly_words_per_exp(Abits, ctx->minfo);
@ -488,6 +540,121 @@ void fmpz_mpoly_from_univar(fmpz_mpoly_t A, const fmpz_mpoly_univar_t B,
}
TMP_END;
fmpz_mpoly_from_univar_bits(A, bits, B, var, ctx);
_fmpz_mpoly_from_univar(A, bits, B, var, ctx);
}
#define COEFF(A, i) ((void*)(A->coeffs + (i)*R->elem_size))
static void mpoly_univar_set_fmpz_mpoly_univar(
mpoly_univar_t A,
mpoly_void_ring_t R,
const fmpz_mpoly_univar_t B,
const fmpz_mpoly_ctx_t ctx)
{
slong i;
mpoly_univar_fit_length(A, B->length, R);
A->length = B->length;
for (i = B->length - 1; i >= 0; i--)
{
fmpz_set(A->exps + i, B->exps + i);
fmpz_mpoly_set(COEFF(A, i), B->coeffs + i, ctx);
}
}
static void mpoly_univar_swap_fmpz_mpoly_univar(
mpoly_univar_t A,
mpoly_void_ring_t R,
fmpz_mpoly_univar_t B,
const fmpz_mpoly_ctx_t ctx)
{
slong i;
mpoly_univar_fit_length(A, B->length, R);
fmpz_mpoly_univar_fit_length(B, A->length, ctx);
for (i = FLINT_MAX(A->length, B->length) - 1; i >= 0; i--)
{
fmpz_swap(A->exps + i, B->exps + i);
fmpz_mpoly_swap(COEFF(A, i), B->coeffs + i, ctx);
}
SLONG_SWAP(A->length, B->length);
}
int fmpz_mpoly_univar_pseudo_gcd(
fmpz_mpoly_univar_t gx,
const fmpz_mpoly_univar_t ax,
const fmpz_mpoly_univar_t bx,
const fmpz_mpoly_ctx_t ctx)
{
int success;
mpoly_void_ring_t R;
mpoly_univar_t Ax, Bx, Gx;
mpoly_void_ring_init_fmpz_mpoly_ctx(R, ctx);
mpoly_univar_init(Ax, R);
mpoly_univar_init(Bx, R);
mpoly_univar_init(Gx, R);
mpoly_univar_set_fmpz_mpoly_univar(Ax, R, ax, ctx);
mpoly_univar_set_fmpz_mpoly_univar(Bx, R, bx, ctx);
success = mpoly_univar_pseudo_gcd_ducos(Gx, Ax, Bx, R);
if (success)
mpoly_univar_swap_fmpz_mpoly_univar(Gx, R, gx, ctx);
mpoly_univar_clear(Ax, R);
mpoly_univar_clear(Bx, R);
mpoly_univar_clear(Gx, R);
return success;
}
int fmpz_mpoly_univar_resultant(
fmpz_mpoly_t d,
const fmpz_mpoly_univar_t ax,
const fmpz_mpoly_univar_t bx,
const fmpz_mpoly_ctx_t ctx)
{
int success;
mpoly_void_ring_t R;
mpoly_univar_t Ax, Bx;
mpoly_void_ring_init_fmpz_mpoly_ctx(R, ctx);
mpoly_univar_init(Ax, R);
mpoly_univar_init(Bx, R);
mpoly_univar_set_fmpz_mpoly_univar(Ax, R, ax, ctx);
mpoly_univar_set_fmpz_mpoly_univar(Bx, R, bx, ctx);
success = mpoly_univar_resultant(d, Ax, Bx, R);
mpoly_univar_clear(Ax, R);
mpoly_univar_clear(Bx, R);
return success;
}
int fmpz_mpoly_univar_discriminant(
fmpz_mpoly_t d,
const fmpz_mpoly_univar_t fx,
const fmpz_mpoly_ctx_t ctx)
{
int success;
mpoly_void_ring_t R;
mpoly_univar_t Fx;
mpoly_void_ring_init_fmpz_mpoly_ctx(R, ctx);
mpoly_univar_init(Fx, R);
mpoly_univar_set_fmpz_mpoly_univar(Fx, R, fx, ctx);
success = mpoly_univar_discriminant(d, Fx, R);
mpoly_univar_clear(Fx, R);
return success;
}

21
fmpz_mpoly_factor.h
View file

@ -35,21 +35,6 @@
/*****************************************************************************/
FLINT_DLL void fmpz_mpoly_fit_length_set_bits(fmpz_mpoly_t A, slong len,
flint_bitcnt_t bits, const fmpz_mpoly_ctx_t ctx);
FLINT_DLL int fmpz_mpoly_is_fmpz_poly(const fmpz_mpoly_t A, slong var,
const fmpz_mpoly_ctx_t ctx);
FLINT_DLL int fmpz_mpoly_get_fmpz_poly(fmpz_poly_t A, const fmpz_mpoly_t B,
slong var, const fmpz_mpoly_ctx_t ctx);
FLINT_DLL void _fmpz_mpoly_set_fmpz_poly(fmpz_mpoly_t A, flint_bitcnt_t Abits,
const fmpz * Bcoeffs, slong Blen, slong var, const fmpz_mpoly_ctx_t ctx);
FLINT_DLL void fmpz_mpoly_set_fmpz_poly(fmpz_mpoly_t A, const fmpz_poly_t B,
slong v, const fmpz_mpoly_ctx_t ctx);
FLINT_DLL void tuple_print(fmpz * alpha, slong n);
FLINT_DLL void tuple_saturate(fmpz * alpha, slong n, slong m);
@ -306,6 +291,12 @@ FLINT_DLL void fmpz_mpoly_from_mpolyv(fmpz_mpoly_t A, flint_bitcnt_t Abits,
FLINT_DLL int _fmpz_mpoly_vec_content_mpoly(fmpz_mpoly_t g,
const fmpz_mpoly_struct * A, slong Alen, const fmpz_mpoly_ctx_t ctx);
FLINT_DLL void _fmpz_mpoly_vec_divexact_mpoly(fmpz_mpoly_struct * A,
slong Alen, const fmpz_mpoly_t c, const fmpz_mpoly_ctx_t ctx);
FLINT_DLL void _fmpz_mpoly_vec_mul_mpoly(fmpz_mpoly_struct * A,
slong Alen, const fmpz_mpoly_t c, const fmpz_mpoly_ctx_t ctx);
/*****************************************************************************/
FLINT_DLL int _fmpz_mpoly_gcd_algo(fmpz_mpoly_t G, fmpz_mpoly_t Abar,

2
fmpz_mpoly_factor/factor.c
View file

@ -394,7 +394,7 @@ static int _refine_content_factors(
{
fmpz_mpolyv_fit_length(g, g->length + 1, ctx);
c = g->coeffs + g->length;
fmpz_mpoly_from_univar_bits(c, bits, u, v, ctx);
_fmpz_mpoly_from_univar(c, bits, u, v, ctx);
g->length++;
}
else

2
fmpz_mpoly_factor/lcc_kaltofen.c
View file

@ -29,7 +29,7 @@ static void fmpz_mpoly_convert_perm(
NA = mpoly_words_per_exp(Abits, Actx->minfo);
NB = mpoly_words_per_exp(B->bits, Bctx->minfo);
fmpz_mpoly_fit_length_set_bits(A, B->length, Abits, Actx);
fmpz_mpoly_fit_length_reset_bits(A, B->length, Abits, Actx);
A->length = B->length;
for (i = 0; i < B->length; i++)
{

23
fmpz_mpoly_factor/mpolyv.c
View file

@ -189,3 +189,26 @@ int _fmpz_mpoly_vec_content_mpoly(
return 1;
}
void _fmpz_mpoly_vec_divexact_mpoly(
fmpz_mpoly_struct * A, slong Alen,
const fmpz_mpoly_t c,
const fmpz_mpoly_ctx_t ctx)
{
slong i;
for (i = 0; i < Alen; i++)
fmpz_mpoly_divexact(A + i, A + i, c, ctx);
}
void _fmpz_mpoly_vec_mul_mpoly(
fmpz_mpoly_struct * A, slong Alen,
const fmpz_mpoly_t c,
const fmpz_mpoly_ctx_t ctx)
{
slong i;
for (i = 0; i < Alen; i++)
fmpz_mpoly_mul(A + i, A + i, c, ctx);
}

247
fmpz_mpoly_factor/test/t-gcd_subresultant.c Normal file
View file

@ -0,0 +1,247 @@
/*
Copyright (C) 2021 Daniel Schultz
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 2.1 of the License, or
(at your option) any later version. See <https://www.gnu.org/licenses/>.
*/
#include "fmpz_mpoly_factor.h"
int compute_gcd(
fmpz_mpoly_t G,
const fmpz_mpoly_t A,
const fmpz_mpoly_t B,
const fmpz_mpoly_ctx_t ctx)
{
int success;
slong var = 0;
fmpz_mpoly_t Ac, Bc, Gc, s, t;
fmpz_mpoly_univar_t Ax, Bx, Gx;
if (fmpz_mpoly_is_zero(A, ctx) || fmpz_mpoly_is_zero(B, ctx))
return fmpz_mpoly_gcd(G, A, B, ctx);
fmpz_mpoly_init(Ac, ctx);
fmpz_mpoly_init(Bc, ctx);
fmpz_mpoly_init(Gc, ctx);
fmpz_mpoly_init(s, ctx);
fmpz_mpoly_init(t, ctx);
fmpz_mpoly_univar_init(Ax, ctx);
fmpz_mpoly_univar_init(Bx, ctx);
fmpz_mpoly_univar_init(Gx, ctx);
fmpz_mpoly_to_univar(Ax, A, var, ctx);
fmpz_mpoly_to_univar(Bx, B, var, ctx);
success = _fmpz_mpoly_vec_content_mpoly(Ac, Ax->coeffs, Ax->length, ctx) &&
_fmpz_mpoly_vec_content_mpoly(Bc, Bx->coeffs, Bx->length, ctx) &&
fmpz_mpoly_gcd(Gc, Ac, Bc, ctx);
if (!success)
goto cleanup;
_fmpz_mpoly_vec_divexact_mpoly(Ax->coeffs, Ax->length, Ac, ctx);
_fmpz_mpoly_vec_divexact_mpoly(Bx->coeffs, Bx->length, Bc, ctx);
success = fmpz_cmp(Ax->exps + 0, Bx->exps + 0) > 0 ?
fmpz_mpoly_univar_pseudo_gcd(Gx, Ax, Bx, ctx) :
fmpz_mpoly_univar_pseudo_gcd(Gx, Bx, Ax, ctx);
if (!success)
goto cleanup;
if (fmpz_mpoly_gcd(t, Ax->coeffs + 0, Bx->coeffs + 0, ctx) &&
t->length == 1)
{
fmpz_mpoly_term_content(s, Gx->coeffs + 0, ctx);
fmpz_mpoly_divexact(t, Gx->coeffs + 0, s, ctx);
_fmpz_mpoly_vec_divexact_mpoly(Gx->coeffs, Gx->length, t, ctx);
}
else if (fmpz_mpoly_gcd(t, Ax->coeffs + Ax->length - 1,
Bx->coeffs + Bx->length - 1, ctx) && t->length == 1)
{
fmpz_mpoly_term_content(s, Gx->coeffs + Gx->length - 1, ctx);
fmpz_mpoly_divexact(t, Gx->coeffs + Gx->length - 1, s, ctx);
_fmpz_mpoly_vec_divexact_mpoly(Gx->coeffs, Gx->length, t, ctx);
}
success = _fmpz_mpoly_vec_content_mpoly(t, Gx->coeffs, Gx->length, ctx);
if (!success)
goto cleanup;
_fmpz_mpoly_vec_divexact_mpoly(Gx->coeffs, Gx->length, t, ctx);
_fmpz_mpoly_vec_mul_mpoly(Gx->coeffs, Gx->length, Gc, ctx);
_fmpz_mpoly_from_univar(G, FLINT_MIN(A->bits, B->bits), Gx, var, ctx);
fmpz_mpoly_unit_normalize(G, ctx);
success = 1;
cleanup:
fmpz_mpoly_clear(Ac, ctx);
fmpz_mpoly_clear(Bc, ctx);
fmpz_mpoly_clear(Gc, ctx);
fmpz_mpoly_clear(s, ctx);
fmpz_mpoly_clear(t, ctx);
fmpz_mpoly_univar_clear(Ax, ctx);
fmpz_mpoly_univar_clear(Bx, ctx);
fmpz_mpoly_univar_clear(Gx, ctx);
return success;
}
void gcd_check(
fmpz_mpoly_t g,
fmpz_mpoly_t a,
fmpz_mpoly_t b,
const fmpz_mpoly_t gdiv,
fmpz_mpoly_ctx_t ctx,
slong i,
slong j,
const char * name)
{
int res;
fmpz_mpoly_t ca, cb, cg;
fmpz_mpoly_init(ca, ctx);
fmpz_mpoly_init(cb, ctx);
fmpz_mpoly_init(cg, ctx);
res = compute_gcd(g, a, b, ctx);
fmpz_mpoly_assert_canonical(g, ctx);
if (!res)
{
flint_printf("FAIL: Check gcd can be computed\n");
flint_printf("i = %wd, j = %wd, %s\n", i, j, name);
flint_abort();
}
if (!fmpz_mpoly_is_zero(gdiv, ctx))
{
if (!fmpz_mpoly_divides(ca, g, gdiv, ctx))
{
flint_printf("FAIL: Check divisor of gcd\n");
flint_printf("i = %wd, j = %wd, %s\n", i, j, name);
flint_abort();
}
}
if (fmpz_mpoly_is_zero(g, ctx))
{
if (!fmpz_mpoly_is_zero(a, ctx) || !fmpz_mpoly_is_zero(b, ctx))
{
flint_printf("FAIL: Check zero gcd\n");
flint_printf("i = %wd, j = %wd, %s\n", i, j, name);
flint_abort();
}
goto cleanup;
}
if (fmpz_sgn(g->coeffs + 0) <= 0)
{
flint_printf("FAIL: Check gcd has positive lc\n");
flint_printf("i = %wd, j = %wd, %s\n", i, j, name);
flint_abort();
}
res = 1;
res = res && fmpz_mpoly_divides(ca, a, g, ctx);
res = res && fmpz_mpoly_divides(cb, b, g, ctx);
if (!res)
{
flint_printf("FAIL: Check divisibility\n");
flint_printf("i = %wd, j = %wd, %s\n", i, j, name);
flint_abort();
}
res = compute_gcd(cg, ca, cb, ctx);
fmpz_mpoly_assert_canonical(cg, ctx);
if (!res)
{
flint_printf("FAIL: Check gcd of cofactors can be computed\n");
flint_printf("i = %wd, j = %wd, %s\n", i, j, name);
flint_abort();
}
if (!fmpz_mpoly_is_one(cg, ctx))
{
flint_printf("FAIL: Check gcd of cofactors is one\n");
flint_printf("i = %wd, j = %wd, %s\n", i, j, name);
flint_abort();
}
cleanup:
fmpz_mpoly_clear(ca, ctx);
fmpz_mpoly_clear(cb, ctx);
fmpz_mpoly_clear(cg, ctx);
}
int
main(void)
{
slong i, j, tmul = 10;
FLINT_TEST_INIT(state);
flint_printf("gcd_subresultant....");
fflush(stdout);
for (i = 0; i < tmul * flint_test_multiplier(); i++)
{
fmpz_mpoly_ctx_t ctx;
fmpz_mpoly_t a, b, g, t;
slong len, len1, len2;
slong degbound;
flint_bitcnt_t coeff_bits;
fmpz_mpoly_ctx_init_rand(ctx, state, 3);
fmpz_mpoly_init(g, ctx);
fmpz_mpoly_init(a, ctx);
fmpz_mpoly_init(b, ctx);
fmpz_mpoly_init(t, ctx);
len = n_randint(state, 15) + 1;
len1 = n_randint(state, 15);
len2 = n_randint(state, 15);
degbound = 1 + 10/ctx->minfo->nvars;
for (j = 0; j < 4; j++)
{
coeff_bits = n_randint(state, 100) + 1;
fmpz_mpoly_randtest_bound(t, state, len, coeff_bits, degbound, ctx);
if (fmpz_mpoly_is_zero(t, ctx))
fmpz_mpoly_one(t, ctx);
coeff_bits = n_randint(state, 100) + 1;
fmpz_mpoly_randtest_bound(a, state, len1, coeff_bits, degbound, ctx);
coeff_bits = n_randint(state, 100) + 1;
fmpz_mpoly_randtest_bound(b, state, len2, coeff_bits, degbound, ctx);
fmpz_mpoly_mul(a, a, t, ctx);
fmpz_mpoly_mul(b, b, t, ctx);
fmpz_mpoly_randtest_bits(g, state, len, coeff_bits, FLINT_BITS, ctx);
gcd_check(g, a, b, t, ctx, i, j, "random dense");
}
fmpz_mpoly_clear(g, ctx);
fmpz_mpoly_clear(a, ctx);
fmpz_mpoly_clear(b, ctx);
fmpz_mpoly_clear(t, ctx);
fmpz_mpoly_ctx_clear(ctx);
}
flint_printf("PASS\n");
FLINT_TEST_CLEANUP(state);
return 0;
}

10
mpoly.h
View file

@ -1090,6 +1090,16 @@ FLINT_DLL flint_bitcnt_t mpoly_exp_bits_required_ffmpz(const fmpz * user_exp,
FLINT_DLL flint_bitcnt_t mpoly_exp_bits_required_pfmpz(fmpz * const * user_exp,
const mpoly_ctx_t mctx);
MPOLY_INLINE
flint_bitcnt_t mpoly_gen_pow_exp_bits_required(slong v, ulong e,
const mpoly_ctx_t mctx)
{
return 1 + FLINT_BIT_COUNT(e); /* only lex and deg supported */
}
FLINT_DLL int mpoly_is_poly(const ulong * Aexps, slong Alen,
flint_bitcnt_t Abits, slong var, const mpoly_ctx_t mctx);
FLINT_DLL void mpoly_pack_vec_ui(ulong * exp1, const ulong * exp2,
flint_bitcnt_t bits, slong nfields, slong len);
FLINT_DLL void mpoly_pack_vec_fmpz(ulong * exp1, const fmpz * exp2,

55
mpoly/is_poly.c Normal file
View file

@ -0,0 +1,55 @@
/*
Copyright (C) 2018-2020 Daniel Schultz
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 2.1 of the License, or
(at your option) any later version. See <https://www.gnu.org/licenses/>.
*/
#include "mpoly.h"
/* TODO: this could a bit faster */
int mpoly_is_poly(
const ulong * Aexps,
slong Alen,
flint_bitcnt_t Abits,
slong var,
const mpoly_ctx_t mctx)
{
int ret = 1;
slong i, j;
slong N = mpoly_words_per_exp(Abits, mctx);
slong nvars = mctx->nvars;
fmpz * t;
TMP_INIT;
TMP_START;
t = (fmpz *) TMP_ALLOC(nvars*sizeof(fmpz));
for (i = 0; i < nvars; i++)
fmpz_init(t + i);
for (i = 0; i < Alen; i++)
{
mpoly_get_monomial_ffmpz(t, Aexps + N*i, Abits, mctx);
for (j = 0; j < nvars; j++)
{
if (j != var && !fmpz_is_zero(t + j))
{
ret = 0;
goto cleanup;
}
}
}
cleanup:
for (i = 0; i < nvars; i++)
fmpz_clear(t + i);
TMP_END;
return ret;
}

16
nmod_mpoly_factor.h
View file

@ -58,13 +58,6 @@ FLINT_DLL int nmod_mat_is_reduced(const nmod_mat_t N);
FLINT_DLL void nmod_mat_init_nullspace_tr(nmod_mat_t X, nmod_mat_t tmp);
FLINT_DLL int mpoly_is_poly(
const ulong * Aexps,
slong Alen,
flint_bitcnt_t Abits,
slong var,
const mpoly_ctx_t mctx);
FLINT_DLL int nmod_mpoly_is_nmod_poly(
const nmod_mpoly_t A,
slong var,
@ -383,15 +376,6 @@ FLINT_DLL void nmod_mpoly_compression_undo(nmod_mpoly_t A,
/*****************************************************************************/
NMOD_MPOLY_FACTOR_INLINE
flint_bitcnt_t mpoly_gen_pow_exp_bits_required(slong v, ulong e,
const mpoly_ctx_t mctx)
{
return 1 + FLINT_BIT_COUNT(e); /* only lex and deg supported */
}
/*****************************************************************************/
FLINT_DLL int nmod_mpolyu_is_canonical(const nmod_mpolyu_t A,
const nmod_mpoly_ctx_t ctx);

44
nmod_mpoly_factor/get_set_is_nmod_poly.c
View file

@ -11,50 +11,6 @@
#include "nmod_mpoly_factor.h"
/* TODO: this could a bit faster */
int mpoly_is_poly(
const ulong * Aexps,
slong Alen,
flint_bitcnt_t Abits,
slong var,
const mpoly_ctx_t mctx)
{
int ret = 1;
slong i, j;
slong N = mpoly_words_per_exp(Abits, mctx);
slong nvars = mctx->nvars;
fmpz * t;
TMP_INIT;
TMP_START;
t = (fmpz *) TMP_ALLOC(nvars*sizeof(fmpz));
for (i = 0; i < nvars; i++)
fmpz_init(t + i);
for (i = 0; i < Alen; i++)
{
mpoly_get_monomial_ffmpz(t, Aexps + N*i, Abits, mctx);
for (j = 0; j < nvars; j++)
{
if (j != var && !fmpz_is_zero(t + j))
{
ret = 0;
goto cleanup;
}
}
}
cleanup:
for (i = 0; i < nvars; i++)
fmpz_clear(t + i);
TMP_END;
return ret;
}
int nmod_mpoly_is_nmod_poly(
const nmod_mpoly_t A,
slong var,