flint3-sys 3.6.0

/*
    Copyright (C) 2025 Marc Mezzarobba

    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 <https://www.gnu.org/licenses/>.
*/

/* todo: make it possible to compute a superset of the result, implement the
 * algorithm of Gerhard, Giesbrecht, Storjohann, Zima (2003) */

#include "fmpz.h"
#include "fmpz_vec.h"
#include "gr.h"
#include "gr_poly.h"
#include "gr_vec.h"

int
gr_poly_dispersion_resultant(fmpz_t disp, gr_vec_t disp_set,
			     const gr_poly_t f, const gr_poly_t g,
			     gr_ctx_t ctx)
{
    gr_ctx_t Pol_a, ZZ;
    gr_poly_t a, fa, ga, res;
    gr_vec_t roots;
    fmpz_vec_t mult;

    truth_t is_zero = truth_and(gr_poly_is_zero(f, ctx),
				gr_poly_is_zero(g, ctx));
    if (is_zero == T_TRUE)
        return GR_DOMAIN;
    else if (is_zero == T_UNKNOWN)
        return GR_UNABLE;

    gr_ctx_init_gr_poly(Pol_a, ctx);
    gr_ctx_init_fmpz(ZZ);

    gr_poly_init(a, ctx);
    gr_poly_init(fa, Pol_a);
    gr_poly_init(ga, Pol_a);
    gr_poly_init(res, ZZ);

    gr_vec_init(roots, 0, ZZ);
    fmpz_vec_init(mult, 0);

    int status = GR_SUCCESS;

    status |= gr_poly_set_gr_poly_other(fa, f, ctx, Pol_a);
    status |= gr_poly_set_gr_poly_other(ga, g, ctx, Pol_a);

    status |= gr_poly_gen(a, ctx);
    status |= gr_poly_taylor_shift(fa, fa, a, Pol_a);

    /* todo: composed sum */
    status |= gr_poly_resultant(res, fa, ga, Pol_a);

    status |= gr_poly_roots_other(roots, mult, res, ctx, 0, ZZ);

    if (disp_set != NULL)
        gr_vec_set_length(disp_set, 0, ZZ);

    for (slong i = 0; i < roots->length; i++)
    {
        gr_srcptr rt = gr_vec_entry_srcptr(roots, i, ZZ);
        if (fmpz_sgn(rt) >= 0)
        {
            if (disp_set != NULL)
                status |= gr_vec_append(disp_set, rt, ZZ);
            if (disp != NULL && fmpz_cmp(rt, disp) >= 0)
                fmpz_swap(disp, (fmpz *) rt);
        }
    }

    if (disp_set != NULL)
        _fmpz_vec_sort(disp_set->entries, disp_set->length);

    fmpz_vec_clear(mult);
    gr_vec_clear(roots, ZZ);

    gr_poly_clear(res, ctx);
    gr_poly_clear(ga, Pol_a);
    gr_poly_clear(fa, Pol_a);
    gr_poly_clear(a, ctx);

    gr_ctx_clear(ZZ);
    gr_ctx_clear(Pol_a);

    return status;
}

int
gr_poly_dispersion_from_factors(fmpz_t disp, gr_vec_t disp_set,
			   const gr_poly_vec_t ffac, const gr_poly_vec_t gfac,
			   gr_ctx_t ctx)
{
    gr_ctx_t ZZ, pctx;
    gr_ptr a, b, _alpha;
    fmpz_t alpha;
    gr_poly_t pshift;

    if (gr_ctx_is_unique_factorization_domain(ctx) != T_TRUE)
        return GR_UNABLE;

    /* todo: generalize to positive characteristic? (cf. Gerhard et al.) */
    if (gr_ctx_is_finite_characteristic(ctx) != T_FALSE)
        return GR_UNABLE;

    gr_ctx_init_fmpz(ZZ);

    if (disp_set != NULL)
        gr_vec_set_length(disp_set, 0, ZZ);

    /* We assume that f and g are nonzero, so no factors means constant */
    if (ffac->length == 0 || gfac->length == 0)
        return GR_SUCCESS;

    gr_ctx_init_gr_poly(pctx, ctx);
    GR_TMP_INIT3(a, b, _alpha, ctx);
    fmpz_init(alpha);
    gr_poly_init(pshift, ctx);

    int status = GR_SUCCESS;

    if (ffac == gfac)
    {
        if (disp_set != NULL)
            status |= gr_vec_append(disp_set, alpha, ZZ);
        if (disp != NULL)
            fmpz_zero(disp);
    }

    for (slong i = 0; i < ffac->length; i++)
    {
        const gr_poly_struct *p = gr_poly_vec_entry_srcptr(ffac, i, ctx);

        slong j0 = (ffac == gfac) ? i + 1 : 0;

        for (slong j = j0; j < gfac->length; j++)
        {
            const gr_poly_struct *q = gr_poly_vec_entry_srcptr(gfac, j, ctx);

            int status1 = GR_SUCCESS;

            status1 |= gr_poly_leading_taylor_shift(_alpha, p, q, ctx);
            if (status1 == GR_DOMAIN)
                continue;
            status1 |= gr_get_fmpz(alpha, _alpha, ctx);
            if (status1 == GR_DOMAIN)
                continue;
            status |= status1;

            if (status != GR_SUCCESS)
                goto cleanup;

            if (ffac == gfac)
                fmpz_abs(alpha, alpha);
            else if (fmpz_sgn(alpha) < 0)
                continue;

            if (disp_set != NULL
                ? gr_vec_contains(disp_set, alpha, ZZ) == T_TRUE
                : disp != NULL && fmpz_cmp(alpha, disp) <= 0)
                continue;

            if (p->length > 2)  /* GR_SUCCESS => degree well-defined */
            {
                status |= gr_poly_taylor_shift(pshift, p, _alpha, ctx);
                int eq = gr_poly_equal(pshift, q, ctx);
                if (eq == T_FALSE)
                    continue;
                else if (eq == T_UNKNOWN)
                    status = GR_UNABLE;
            }

            if (status != GR_SUCCESS)
                goto cleanup;

            if (disp_set != NULL)
                status |= gr_vec_append(disp_set, alpha, ZZ);
            if (disp != NULL && fmpz_cmp(alpha, disp) > 0)
                fmpz_swap(disp, alpha);
        }
    }

    if (disp_set != NULL)
        _fmpz_vec_sort(disp_set->entries, disp_set->length);

cleanup:

    gr_poly_clear(pshift, ctx);
    fmpz_clear(alpha);
    GR_TMP_CLEAR3(a, b, _alpha, ctx);
    gr_ctx_clear(pctx);
    gr_ctx_clear(ZZ);

    return status;
}

int
gr_poly_dispersion_factor(fmpz_t disp, gr_vec_t disp_set,
			  const gr_poly_t f, const gr_poly_t g,
			  gr_ctx_t ctx)
{
    gr_ctx_t pctx, ZZ;
    gr_poly_t fsqf, gsqf;
    gr_poly_vec_t ffac, gfac;
    fmpz_vec_t ignored;
    gr_ptr c;

    switch (truth_and(gr_poly_is_zero(f, ctx), gr_poly_is_zero(g, ctx)))
    {
        case T_TRUE:
            return GR_DOMAIN;
        case T_UNKNOWN:
            return GR_UNABLE;
        case T_FALSE:
            ;
    }

    int status = GR_SUCCESS;

    gr_ctx_init_fmpz(ZZ);
    gr_ctx_init_gr_poly(pctx, ctx);
    gr_poly_init(fsqf, ctx);
    gr_poly_init(gsqf, ctx);
    gr_poly_vec_init(ffac, 0, ctx);
    gr_poly_vec_init(gfac, 0, ctx);
    fmpz_vec_init(ignored, 0);
    GR_TMP_INIT(c, pctx);

    if (gr_poly_squarefree_part(fsqf, f, ctx) != GR_SUCCESS)
        status |= gr_poly_set(fsqf, f, ctx);

    status |= gr_factor(c, (gr_vec_struct *) ffac, ignored, f, 0, pctx);

    if (f == g)
    {
        status |= gr_poly_dispersion_from_factors(disp, disp_set, ffac, ffac, ctx);
    }
    else
    {
        if (gr_poly_squarefree_part(gsqf, g, ctx) != GR_SUCCESS)
            status |= gr_poly_set(gsqf, g, ctx);
        status |= gr_factor(c, (gr_vec_struct *) gfac, ignored, g, 0, pctx);
        status |= gr_poly_dispersion_from_factors(disp, disp_set, ffac, gfac, ctx);
    }

    GR_TMP_CLEAR(c, pctx);
    gr_poly_vec_clear(gfac, ctx);
    gr_poly_vec_clear(ffac, ctx);
    gr_poly_clear(fsqf, ctx);
    gr_poly_clear(gsqf, ctx);
    fmpz_vec_clear(ignored);
    gr_ctx_clear(pctx);
    gr_ctx_clear(ZZ);

    return status;
}

int
gr_poly_dispersion(fmpz_t disp, gr_vec_t disp_set,
		   const gr_poly_t f, const gr_poly_t g,
		   gr_ctx_t ctx)
{
    return gr_poly_dispersion_factor(disp, disp_set, f, g, ctx);
}