flintinterface.cc

Go to the documentation of this file.

 
/* #[ License : */
/*
 *   Copyright (C) 1984-2026 J.A.M. Vermaseren
 *   When using this file you are requested to refer to the publication
 *   J.A.M.Vermaseren "New features of FORM" math-ph/0010025
 *   This is considered a matter of courtesy as the development was paid
 *   for by FOM the Dutch physics granting agency and we would like to
 *   be able to track its scientific use to convince FOM of its value
 *   for the community.
 *
 *   This file is part of FORM.
 *
 *   FORM is free software: you can redistribute it and/or modify it under the
 *   terms of the GNU General Public License as published by the Free Software
 *   Foundation, either version 3 of the License, or (at your option) any later
 *   version.
 *
 *   FORM is distributed in the hope that it will be useful, but WITHOUT ANY
 *   WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 *   FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
 *   details.
 *
 *   You should have received a copy of the GNU General Public License along
 *   with FORM.  If not, see <http://www.gnu.org/licenses/>.
 */
/* #] License : */
 
#ifdef HAVE_CONFIG_H
#ifndef CONFIG_H_INCLUDED
#define CONFIG_H_INCLUDED
#include <config.h>
#endif
#endif
 
extern "C" {
#include "form3.h"
}
 
#include "flintinterface.h"
 
/*
    #[ Types
 * Use fixed-size integer types (int32_t, uint32_t, int64_t, uint64_t) except when we refer
 * specifically to FORM term data, when we use WORD. This code has only been tested on 64bit
 * systems where WORDs are int32_t, and some functions (fmpz_get_form, fmpz_set_form) require this
 * to be the case. Enforce this:
 */
static_assert(sizeof(WORD) == sizeof(int32_t), "flint interface expects 32bit WORD");
static_assert(sizeof(UWORD) == sizeof(uint32_t), "flint interface expects 32bit WORD");
static_assert(BITSINWORD == 32, "flint interface expects 32bit WORD");
static_assert(sizeof(ulong) == sizeof(uint64_t), "flint interface expects ulong is uint64_t");
static_assert(sizeof(slong) == sizeof(int64_t), "flint interface expects slong is int64_t");
 
/*
 * Flint functions take arguments or return values which may be "slong" or "ulong" in its
 * documentation, and these are int64_t and uint64_t respectively.
 *
    #] Types
*/
 
/*
 * FLINT's univariate poly has a dense representation. For sufficiently sparse polynomials it is
 * faster to use mpoly instead, which is sparse. For a density <= this threshold we switch, where
 * the density defined as is "number of terms" / "maximum degree".
 */
#define UNIVARIATE_DENSITY_THR 0.02f
 
/*
    #[ flint::cleanup :
*/
void flint::cleanup(void) {
    flint_cleanup();
}
/*
    #] flint::cleanup :
    #[ flint::cleanup_master :
*/
void flint::cleanup_master(void) {
    flint_cleanup_master();
}
/*
    #] flint::cleanup_master :
 
    #[ flint::divmod_mpoly :
*/
WORD* flint::divmod_mpoly(PHEAD const WORD *a, const WORD *b, const bool return_rem,
    const WORD must_fit_term, const var_map_t &var_map, const bool sort_vars) {
 
    flint::mpoly_ctx ctx(var_map.size());
    flint::mpoly pa(ctx.d), pb(ctx.d), denpa(ctx.d), denpb(ctx.d);
 
    flint::from_argument_mpoly(pa.d, denpa.d, a, false, var_map, ctx.d);
    flint::from_argument_mpoly(pb.d, denpb.d, b, false, var_map, ctx.d);
 
    // The input won't have any symbols with negative powers, but there may be rational
    // coefficients. Verify this:
    if ( fmpz_mpoly_is_fmpz(denpa.d, ctx.d) != 1 ) {
        // INTERNAL_ERROR_EXCL_START
        MLOCK(ErrorMessageLock);
        MesPrint("!>flint::divmod_mpoly: error: denpa is non-constant");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
        // INTERNAL_ERROR_EXCL_STOP
    }
    if ( fmpz_mpoly_is_fmpz(denpb.d, ctx.d) != 1 ) {
        // INTERNAL_ERROR_EXCL_START
        MLOCK(ErrorMessageLock);
        MesPrint("!>flint::divmod_mpoly: error: denpb is non-constant");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
        // INTERNAL_ERROR_EXCL_STOP
    }
 
 
    flint::fmpz scale;
    flint::mpoly div(ctx.d), rem(ctx.d);
 
    fmpz_mpoly_quasidivrem(scale.d, div.d, rem.d, pa.d, pb.d, ctx.d);
 
    // The quotient must be multiplied by the denominator of the divisor
    fmpz_mpoly_mul(div.d, div.d, denpb.d, ctx.d);
 
    // The overall denominator of both div and rem is given by scale*denpa. This we will pass to
    // to_argument_mpoly's "denscale" argument which performs the division in the result. We have
    // already checked denpa is just an fmpz.
    fmpz_mul(scale.d, scale.d, fmpz_mpoly_term_coeff_ref(denpa.d, 0, ctx.d));
 
 
    // Determine the size of the output by passing write = false.
    const bool with_arghead = false;
    bool write = false;
    const uint64_t prev_size = 0;
    const uint64_t out_size = return_rem ?
        (uint64_t)flint::to_argument_mpoly(BHEAD NULL, with_arghead, must_fit_term, write, prev_size,
            rem.d, var_map, ctx.d, sort_vars, scale.d)
        :
        (uint64_t)flint::to_argument_mpoly(BHEAD NULL, with_arghead, must_fit_term, write, prev_size,
            div.d, var_map, ctx.d, sort_vars, scale.d)
        ;
    WORD* res = (WORD *)Malloc1(sizeof(WORD)*out_size, "flint::divrem_mpoly");
 
 
    // Write out the result
    write = true;
    if ( return_rem ) {
        flint::to_argument_mpoly(BHEAD res, with_arghead, must_fit_term, write, prev_size,
            rem.d, var_map, ctx.d, sort_vars, scale.d);
    }
    else {
        flint::to_argument_mpoly(BHEAD res, with_arghead, must_fit_term, write, prev_size,
            div.d, var_map, ctx.d, sort_vars, scale.d);
    }
 
    return res;
}
/*
    #] flint::divmod_mpoly :
    #[ flint::divmod_poly :
*/
WORD* flint::divmod_poly(PHEAD const WORD *a, const WORD *b, const bool return_rem,
    const WORD must_fit_term, const var_map_t &var_map) {
 
    flint::poly pa, pb, denpa, denpb;
 
    flint::from_argument_poly(pa.d, denpa.d, a, false);
    flint::from_argument_poly(pb.d, denpb.d, b, false);
 
    // The input won't have any symbols with negative powers, but there may be rational
    // coefficients. Verify this:
    if ( fmpz_poly_length(denpa.d) != 1 ) {
        // INTERNAL_ERROR_EXCL_START
        MLOCK(ErrorMessageLock);
        MesPrint("!>flint::divmod_poly: error: denpa is non-constant");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
        // INTERNAL_ERROR_EXCL_STOP
    }
    if ( fmpz_poly_length(denpb.d) != 1 ) {
        // INTERNAL_ERROR_EXCL_START
        MLOCK(ErrorMessageLock);
        MesPrint("!>flint::divmod_poly: error: denpb is non-constant");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
        // INTERNAL_ERROR_EXCL_STOP
    }
 
    flint::fmpz scale;
    uint64_t scale_power = 0;
    flint::poly div, rem;
 
    // Get the leading coefficient of pb. fmpz_poly_pseudo_divrem returns only the scaling power.
    fmpz_poly_get_coeff_fmpz(scale.d, pb.d, fmpz_poly_degree(pb.d));
 
    fmpz_poly_pseudo_divrem(div.d, rem.d, (ulong*)&scale_power, pa.d, pb.d);
 
    // The quotient must be multiplied by the denominator of the divisor
    fmpz_poly_mul(div.d, div.d, denpb.d);
 
    // The overall denominator of both div and rem is given by scale^scale_power*denpa. This we will
    // pass to to_argument_poly's "denscale" argument which performs the division in the result. We
    // have already checked denpa is just an fmpz.
    fmpz_pow_ui(scale.d, scale.d, scale_power);
    fmpz_mul(scale.d, scale.d, fmpz_poly_get_coeff_ptr(denpa.d, 0));
 
 
    // Determine the size of the output by passing write = false.
    const bool with_arghead = false;
    bool write = false;
    const uint64_t prev_size = 0;
    const uint64_t out_size = return_rem ?
        (uint64_t)flint::to_argument_poly(BHEAD NULL, with_arghead, must_fit_term, write, prev_size,
            rem.d, var_map, scale.d)
        :
        (uint64_t)flint::to_argument_poly(BHEAD NULL, with_arghead, must_fit_term, write, prev_size,
            div.d, var_map, scale.d)
        ;
    WORD* res = (WORD *)Malloc1(sizeof(WORD)*out_size, "flint::divrem_poly");
 
 
    // Write out the result
    write = true;
    if ( return_rem ) {
        flint::to_argument_poly(BHEAD res, with_arghead, must_fit_term, write, prev_size,
            rem.d, var_map, scale.d);
    }
    else {
        flint::to_argument_poly(BHEAD res, with_arghead, must_fit_term, write, prev_size,
            div.d, var_map, scale.d);
    }
 
    return res;
}
/*
    #] flint::divmod_poly :
 
    #[ flint::factorize_mpoly :
*/
WORD* flint::factorize_mpoly(PHEAD const WORD *argin, WORD *argout, const bool with_arghead,
    const bool is_fun_arg, const var_map_t &var_map, const bool sort_vars) {
 
    flint::mpoly_ctx ctx(var_map.size());
    flint::mpoly arg(ctx.d), den(ctx.d), base(ctx.d);
 
    flint::from_argument_mpoly(arg.d, den.d, argin, with_arghead, var_map, ctx.d);
    // The denominator must be 1:
    if ( fmpz_mpoly_is_one(den.d, ctx.d) != 1 ) {
        // INTERNAL_ERROR_EXCL_START
        MLOCK(ErrorMessageLock);
        MesPrint("!>flint::factorize_mpoly error: den != 1");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
        // INTERNAL_ERROR_EXCL_STOP
    }
 
 
    // Now we can factor the mpoly:
    flint::mpoly_factor arg_fac(ctx.d);
    fmpz_mpoly_factor(arg_fac.d, arg.d, ctx.d);
    const int64_t num_factors = fmpz_mpoly_factor_length(arg_fac.d, ctx.d);
 
    flint::fmpz overall_constant;
    fmpz_mpoly_factor_get_constant_fmpz(overall_constant.d, arg_fac.d, ctx.d);
    // FORM should always have taken the overall constant out in the content. Thus this overall
    // constant factor should be +- 1 here. Verify this:
    if ( ! ( fmpz_equal_si(overall_constant.d, 1) || fmpz_equal_si(overall_constant.d, -1) ) ) {
        // INTERNAL_ERROR_EXCL_START
        MLOCK(ErrorMessageLock);
        MesPrint("!>flint::factorize_mpoly error: overall constant factor != +-1");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
        // INTERNAL_ERROR_EXCL_STOP
    }
 
    // Construct the output. If argout is not NULL, we write the result there.
    // Otherwise, allocate memory.
    // The output is zero-terminated list of factors. If with_arghead, each has an arghead which
    // contains its size. Otherwise, the factors are zero separated.
 
    // We only need to determine the size of the output if we are allocating memory, but we need to
    // loop through the factors to fix their signs anyway. Do both together in one loop:
 
    // Initially 1, for the final trailing 0.
    uint64_t output_size = 1;
 
    // For finding the highest symbol, in FORM's lexicographic ordering
    vector<uint32_t> var_map_inv;
    var_map_inv.resize(var_map.size());
    for ( auto x : var_map ) {
        var_map_inv[x.second] = x.first;
    }
 
    // Store whether we should flip the factor sign in the ouput:
    vector<int32_t> base_sign(num_factors, 0);
 
    for ( int64_t i = 0; i < num_factors; i++ ) {
        fmpz_mpoly_factor_get_base(base.d, arg_fac.d, i, ctx.d);
        const int64_t exponent = fmpz_mpoly_factor_get_exp_si(arg_fac.d, i, ctx.d);
 
        // poly_factorize makes sure the highest power of the "highest symbol" (in FORM's
        // lexicographic ordering) has a positive coefficient. Check this, update overall_constant
        // of the factorization if necessary.
        // Store the sign per factor, so that we can flip the signs in the output without re-checking
        // the individual terms again.
        uint32_t max_var = 0; // FORM symbols start at 20, 0 is a good initial value.
        int32_t max_pow = -1;
        vector<int64_t> base_term_exponents(var_map.size(), 0);
 
        for ( int64_t j = 0; j < fmpz_mpoly_length(base.d, ctx.d); j++ ) {
            fmpz_mpoly_get_term_exp_si((slong*)base_term_exponents.data(), base.d, (slong)j, ctx.d);
 
            for ( size_t k = 0; k < var_map.size(); k++ ) {
                if ( base_term_exponents[k] > 0 && ( var_map_inv[k] > max_var ||
                    ( var_map_inv[k] == max_var && base_term_exponents[k] > max_pow ) ) ) {
 
                    max_var = var_map_inv[k];
                    max_pow = base_term_exponents[k];
                    base_sign[i] = fmpz_sgn(fmpz_mpoly_term_coeff_ref(base.d, j, ctx.d));
                }
            }
        }
        // If this base's sign will be flipped an odd number of times, there is a contribution to
        // the overall sign of the whole factorization:
        if ( ( base_sign[i] == -1 ) && ( exponent % 2 == 1 ) ) {
            fmpz_neg(overall_constant.d, overall_constant.d);
        }
 
        // Now determine the output size of the factor, if we are allocating the memory
        if ( argout == NULL ) {
            const bool write = false;
            for ( int64_t j = 0; j < exponent; j++ ) {
                output_size += (uint64_t)flint::to_argument_mpoly(BHEAD NULL, with_arghead,
                    is_fun_arg, write, 0, base.d, var_map, ctx.d, sort_vars);
            }
        }
    }
    if ( fmpz_sgn(overall_constant.d) == -1 && argout == NULL ) {
        // Add space for a fast-notation number or a normal-notation number and zero separator
        output_size += with_arghead ? 2 : 4+1;
    }
 
    // Now make the allocation if necessary:
    if ( argout == NULL ) {
        argout = (WORD*)Malloc1(sizeof(WORD)*output_size, "flint::factorize_mpoly");
    }
 
 
    // And now comes the actual output:
    WORD* old_argout = argout;
 
    // If the overall sign is negative, first write a full-notation -1. It will be absorbed into the
    // overall factor in the content by the caller.
    if ( fmpz_sgn(overall_constant.d) == -1 ) {
        if ( with_arghead ) {
            // poly writes in fast notation in this case. Fast notation is expected by the caller, to
            // properly merge it with the overall factor of the content.
            *argout++ = -SNUMBER;
            *argout++ = -1;
        }
        else {
            *argout++ = 4; // term size
            *argout++ = 1; // numerator
            *argout++ = 1; // denominator
            *argout++ = -3; // coeff size, negative number
            *argout++ = 0; // factor separator
        }
    }
 
    for ( int64_t i = 0; i < num_factors; i++ ) {
        fmpz_mpoly_factor_get_base(base.d, arg_fac.d, i, ctx.d);
        const int64_t exponent = fmpz_mpoly_factor_get_exp_si(arg_fac.d, i, ctx.d);
 
        if ( base_sign[i] == -1 ) {
            fmpz_mpoly_neg(base.d, base.d, ctx.d);
        }
 
        const bool write = true;
        for ( int64_t j = 0; j < exponent; j++ ) {
            argout += flint::to_argument_mpoly(BHEAD argout, with_arghead, is_fun_arg, write,
                argout-old_argout, base.d, var_map, ctx.d, sort_vars);
        }
    }
    // Final trailing zero to denote the end of the factors.
    *argout++ = 0;
 
    return old_argout;
}
/*
    #] flint::factorize_mpoly :
    #[ flint::factorize_poly :
*/
WORD* flint::factorize_poly(PHEAD const WORD *argin, WORD *argout, const bool with_arghead,
    const bool is_fun_arg, const var_map_t &var_map) {
 
    flint::poly arg, den;
 
    flint::from_argument_poly(arg.d, den.d, argin, with_arghead);
    // The denominator must be 1:
    if ( fmpz_poly_is_one(den.d) != 1 ) {
        // INTERNAL_ERROR_EXCL_START
        MLOCK(ErrorMessageLock);
        MesPrint("!>flint::factorize_poly error: den != 1");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
        // INTERNAL_ERROR_EXCL_STOP
    }
 
 
    // Now we can factor the poly:
    flint::poly_factor arg_fac;
    fmpz_poly_factor(arg_fac.d, arg.d);
    // fmpz_poly_factor_t lacks some convenience functions which fmpz_mpoly_factor_t has.
    // I have worked out how to get the factors by looking at how fmpz_poly_factor_print works.
    const long num_factors = (arg_fac.d)->num;
 
 
    // Construct the output. If argout is not NULL, we write the result there.
    // Otherwise, allocate memory.
    // The output is zero-terminated list of factors. If with_arghead, each has an arghead which
    // contains its size. Otherwise, the factors are zero separated.
    if ( argout == NULL ) {
        // First we need to determine the size of the output. This is the same procedure as the
        // loop below, but we don't write anything in to_argument_poly (write arg: false).
        // Initially 1, for the final trailing 0.
        uint64_t output_size = 1;
 
        // If the overall constant is not 1, make an fmpz_poly out of it and include it. Typically
        // this does not happen: FORM takes out the overall factor. But we do sometimes have this
        // case, for eg awkward interactions with "On highfirst;" and non-dirty arguments.
        if ( ! fmpz_is_one(&(arg_fac.d)->c) ) {
            flint::poly overall_factor;
            fmpz_poly_set_fmpz(overall_factor.d, &(arg_fac.d)->c);
            const bool write = false;
            output_size += (uint64_t)flint::to_argument_poly(BHEAD NULL, with_arghead, is_fun_arg,
                write, 0, overall_factor.d, var_map);
        }
 
        for ( long i = 0; i < num_factors; i++ ) {
            fmpz_poly_struct* base = (arg_fac.d)->p + i;
            
            const long exponent = (arg_fac.d)->exp[i];
 
            const bool write = false;
            for ( long j = 0; j < exponent; j++ ) {
                output_size += (uint64_t)flint::to_argument_poly(BHEAD NULL, with_arghead,
                    is_fun_arg, write, 0, base, var_map);
            }
        }
 
        argout = (WORD*)Malloc1(sizeof(WORD)*output_size, "flint::factorize_poly");
    }
 
    WORD* old_argout = argout;
 
    // If the overall constant is not 1, make an fmpz_poly out of it and include it. Typically
    // this does not happen: FORM takes out the overall factor. But we do sometimes have this
    // case, for eg awkward interactions with "On highfirst;" and non-dirty arguments.
    if ( ! fmpz_is_one(&(arg_fac.d)->c) ) {
        flint::poly overall_factor;
        fmpz_poly_set_fmpz(overall_factor.d, &(arg_fac.d)->c);
        const bool write = true;
        argout += (uint64_t)flint::to_argument_poly(BHEAD argout, with_arghead, is_fun_arg, write,
            argout-old_argout, overall_factor.d, var_map);
    }
 
    for ( long i = 0; i < num_factors; i++ ) {
        fmpz_poly_struct* base = (arg_fac.d)->p + i;
        
        const long exponent = (arg_fac.d)->exp[i];
 
        const bool write = true;
        for ( long j = 0; j < exponent; j++ ) {
            argout += flint::to_argument_poly(BHEAD argout, with_arghead, is_fun_arg, write,
                argout-old_argout, base, var_map);
        }
    }
    *argout = 0;
 
    return old_argout;
}
/*
    #] flint::factorize_poly :
 
    #[ flint::form_sort :
*/
// Sort terms using form's sorting routines. Uses a custom (faster) compare routine, since here
// only symbols can appear.
// This is a modified poly_sort from polywrap.cc.
void flint::form_sort(PHEAD WORD *terms) {
 
    if ( terms[0] < 0 ) {
        // Fast notation, there is nothing to do
        return;
    }
 
    const WORD oldsorttype = AR.SortType;
    AR.SortType = SORTHIGHFIRST;
 
    const WORD in_size = terms[0] - ARGHEAD;
    WORD out_size;
 
    if ( NewSort(BHEAD0) ) {
        Terminate(-1);
    }
    AR.CompareRoutine = (COMPAREDUMMY)(&CompareSymbols);
 
    // Make sure the symbols are in the right order within the terms
    for ( WORD i = ARGHEAD; i < terms[0]; i += terms[i] ) {
        if ( SymbolNormalize(terms+i) < 0 || StoreTerm(BHEAD terms+i) ) {
            AR.SortType = oldsorttype;
            AR.CompareRoutine = (COMPAREDUMMY)(&Compare1);
            LowerSortLevel();
            Terminate(-1);
        }
    }
 
    if ( ( out_size = EndSort(BHEAD terms+ARGHEAD, 1) ) < 0 ) {
        AR.SortType = oldsorttype;
        AR.CompareRoutine = (COMPAREDUMMY)(&Compare1);
        Terminate(-1);
    }
 
    // Check the final size
    if ( in_size != out_size  ) {
        // INTERNAL_ERROR_EXCL_START
        MLOCK(ErrorMessageLock);
        MesPrint("!>flint::form_sort: error: unexpected sorted arg length change %d->%d", in_size,
            out_size);
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
        // INTERNAL_ERROR_EXCL_STOP
    }
 
    AR.SortType = oldsorttype;
    AR.CompareRoutine = (COMPAREDUMMY)(&Compare1);
    terms[1] = 0; // set dirty flag to zero
}
/*
    #] flint::form_sort :
 
    #[ flint::from_argument_mpoly :
*/
// Convert a FORM argument (or 0-terminated list of terms: with_arghead == false) to a
// (multi-variate) fmpz_mpoly_t poly. The "denominator" is return in denpoly, and contains the
// overall negative-power numeric and symbolic factor.
uint64_t flint::from_argument_mpoly(fmpz_mpoly_t poly, fmpz_mpoly_t denpoly, const WORD *args,
    const bool with_arghead, const var_map_t &var_map, const fmpz_mpoly_ctx_t ctx) {
 
    // Some callers re-use their poly, denpoly to avoid calling init/clear unnecessarily.
    // Make sure they are 0 to begin.
    fmpz_mpoly_set_si(poly, 0, ctx);
    fmpz_mpoly_set_si(denpoly, 0, ctx);
 
    // First check for "fast notation" arguments:
    if ( *args == -SNUMBER ) {
        fmpz_mpoly_set_si(poly, *(args+1), ctx);
        fmpz_mpoly_set_si(denpoly, 1, ctx);
        return 2;
    }
 
    if ( *args == -SYMBOL ) {
        // A "fast notation" SYMBOL has a power and coefficient of 1:
        vector<uint64_t> exponents(var_map.size(), 0);
        exponents[var_map.at(*(args+1))] = 1;
 
        fmpz_mpoly_set_coeff_ui_ui(poly, (ulong)1, (ulong*)exponents.data(), ctx);
        fmpz_mpoly_set_ui(denpoly, (ulong)1, ctx);
        return 2;
    }
 
 
    // Now we can iterate through the terms of the argument. If we have
    // an ARGHEAD, we already know where to terminate. Otherwise we'll have
    // to loop until the terminating 0.
    const WORD* arg_stop = with_arghead ? args+args[0] : (WORD*)UINT64_MAX;
    uint64_t arg_size = 0;
    if ( with_arghead ) {
        arg_size = args[0];
        args += ARGHEAD;
    }
 
 
    // Search for numerical or symbol denominators to create "denpoly".
    flint::fmpz den_coeff, tmp;
    fmpz_set_si(den_coeff.d, 1);
    vector<uint64_t> exponents(var_map.size(), 0);
    vector<uint64_t> neg_exponents(var_map.size(), 0);
 
    for ( const WORD* term = args; term < arg_stop; term += term[0] ) {
        const WORD* term_stop = term+term[0];
        const WORD  coeff_size = (term_stop)[-1];
        const WORD* symbol_stop = term_stop - ABS(coeff_size);
        const WORD* t = term;
 
        t++;
        if ( t == symbol_stop ) {
            // Just a number, no symbols
        }
        else {
            t++; // this entry is SYMBOL
            t++; // this entry just has the size of the symbol array, but we can use symbol_stop
 
            for ( const WORD* s = t; s < symbol_stop; s += 2 ) {
                if ( *(s+1) < 0 ) {
                    neg_exponents[var_map.at(*s)] =
                        MaX(neg_exponents[var_map.at(*s)], (uint64_t)(-(*(s+1))) );
                }
            }
        }
 
        // Now check for a denominator in the coefficient:
        if ( *(symbol_stop+ABS(coeff_size/2)) != 1 ) {
            flint::fmpz_set_form(tmp.d, (UWORD*)(symbol_stop+ABS(coeff_size/2)), ABS(coeff_size/2));
            // Record the LCM of the coefficient denominators:
            fmpz_lcm(den_coeff.d, den_coeff.d, tmp.d);
        }
 
        if ( !with_arghead && *term_stop == 0 ) {
            // + 1 for the terminating 0
            arg_size = term_stop - args + 1;
            break;
        }
 
    }
    // Assemble denpoly.
    fmpz_mpoly_set_coeff_fmpz_ui(denpoly, den_coeff.d, (ulong*)neg_exponents.data(), ctx);
 
 
    // For the term coefficients
    flint::fmpz coeff;
 
    for ( const WORD* term = args; term < arg_stop; term += term[0] ) {
        const WORD* term_stop = term+term[0];
        const WORD  coeff_size = (term_stop)[-1];
        const WORD* symbol_stop = term_stop - ABS(coeff_size);
        const WORD* t = term;
 
        fill(exponents.begin(), exponents.end(), 0);
 
        t++; // skip over the total size entry
        if ( t == symbol_stop ) {
            // Just a number, no symbols
        }
        else {
            t++; // this entry is SYMBOL
            t++; // this entry just has the size of the symbol array, but we can use symbol_stop
            for ( const WORD* s = t; s < symbol_stop; s += 2 ) {
                exponents[var_map.at(*s)] = *(s+1);
            }
        }
        // Now read the coefficient
        flint::fmpz_set_form(coeff.d, (UWORD*)symbol_stop, coeff_size/2);
 
        // Multiply by denominator LCM
        fmpz_mul(coeff.d, coeff.d, den_coeff.d);
 
        // Shift by neg_exponents
        for ( size_t i = 0; i < var_map.size(); i++ ) {
            exponents[i] += neg_exponents[i];
        }
 
        // Read the denominator if there is one, and divide it out of the coefficient
        if ( *(symbol_stop+ABS(coeff_size/2)) != 1 ) {
            flint::fmpz_set_form(tmp.d, (UWORD*)(symbol_stop+ABS(coeff_size/2)), ABS(coeff_size/2));
            // By construction, this is an exact division
            fmpz_divexact(coeff.d, coeff.d, tmp.d);
        }
 
        // Push the term to the mpoly, remember to sort when finished! This is much faster than using
        // fmpz_mpoly_set_coeff_fmpz_ui when the terms arrive in the "wrong order".
        fmpz_mpoly_push_term_fmpz_ui(poly, coeff.d, (ulong*)exponents.data(), ctx);
 
        if ( !with_arghead && *term_stop == 0 ) {
            break;
        }
 
    }
    // And now sort the mpoly if necessary:
    if ( ! fmpz_mpoly_is_canonical(poly, ctx) ) {
        fmpz_mpoly_sort_terms(poly, ctx);
    }
 
    return arg_size;
}
/*
    #] flint::from_argument_mpoly :
    #[ flint::from_argument_poly :
*/
// Convert a FORM argument (or 0-terminated list of terms: with_arghead == false) to a
// (uni-variate) fmpz_poly_t poly. The "denominator" is return in denpoly, and contains the
// overall negative-power numeric and symbolic factor.
uint64_t flint::from_argument_poly(fmpz_poly_t poly, fmpz_poly_t denpoly, const WORD *args,
    const bool with_arghead) {
 
    // Some callers re-use their poly, denpoly to avoid calling init/clear unnecessarily.
    // Make sure they are 0 to begin.
    fmpz_poly_set_si(poly, 0);
    fmpz_poly_set_si(denpoly, 0);
 
    // First check for "fast notation" arguments:
    if ( *args == -SNUMBER ) {
        fmpz_poly_set_si(poly, *(args+1));
        fmpz_poly_set_si(denpoly, 1);
        return 2;
    }
    
    if ( *args == -SYMBOL ) {
        // A "fast notation" SYMBOL has a power and coefficient of 1:
        fmpz_poly_set_coeff_si(poly, 1, 1);
        fmpz_poly_set_si(denpoly, 1);
        return 2;
    }
 
    // Now we can iterate through the terms of the argument. If we have
    // an ARGHEAD, we already know where to terminate. Otherwise we'll have
    // to loop until the terminating 0.
    const WORD* arg_stop = with_arghead ? args+args[0] : (WORD*)UINT64_MAX;
    uint64_t arg_size = 0;
    if ( with_arghead ) {
        arg_size = args[0];
        args += ARGHEAD;
    }
 
    // Search for numerical or symbol denominators to create "denpoly".
    flint::fmpz den_coeff, tmp;
    fmpz_set_si(den_coeff.d, 1);
    uint64_t neg_exponent = 0;
 
    for ( const WORD* term = args; term < arg_stop; term += term[0] ) {
 
        const WORD* term_stop = term+term[0];
        const WORD  coeff_size = (term_stop)[-1];
        const WORD* symbol_stop = term_stop - ABS(coeff_size);
        const WORD* t = term;
 
        t++; // skip over the total size entry
        if ( t == symbol_stop ) {
            // Just a number, no symbols
        }
        else {
            t++; // this entry is SYMBOL
            t++; // this entry is the size of the symbol array
            t++; // this is the first (and only) symbol code
            if ( *t < 0 ) {
                neg_exponent = MaX(neg_exponent, (uint64_t)(-(*t)) );
            }
        }
 
        // Now check for a denominator in the coefficient:
        if ( *(symbol_stop+ABS(coeff_size/2)) != 1 ) {
            flint::fmpz_set_form(tmp.d, (UWORD*)(symbol_stop+ABS(coeff_size/2)), ABS(coeff_size/2));
            // Record the LCM of the coefficient denominators:
            fmpz_lcm(den_coeff.d, den_coeff.d, tmp.d);
        }
 
        if ( *term_stop == 0 ) {
            // + 1 for the terminating 0
            arg_size = term_stop - args + 1;
            break;
        }
    }
    // Assemble denpoly.
    fmpz_poly_set_coeff_fmpz(denpoly, neg_exponent, den_coeff.d);
 
 
    // For the term coefficients
    flint::fmpz coeff;
 
    for ( const WORD* term = args; term < arg_stop; term += term[0] ) {
 
        const WORD* term_stop = term+term[0];
        const WORD  coeff_size = (term_stop)[-1];
        const WORD* symbol_stop = term_stop - ABS(coeff_size);
        const WORD* t = term;
 
        uint64_t exponent = 0;
 
        t++; // skip over the total size entry
        if ( t == symbol_stop ) {
            // Just a number, no symbols
        }
        else {
            t++; // this entry is SYMBOL
            t++; // this entry is the size of the symbol array
            t++; // this is the first (and only) symbol code
            exponent = *t++;
        }
 
        // Now read the coefficient
        flint::fmpz_set_form(coeff.d, (UWORD*)symbol_stop, coeff_size/2);
 
        // Multiply by denominator LCM
        fmpz_mul(coeff.d, coeff.d, den_coeff.d);
 
        // Shift by neg_exponent
        exponent += neg_exponent;
 
        // Read the denominator if there is one, and divide it out of the coefficient
        if ( *(symbol_stop+ABS(coeff_size/2)) != 1 ) {
            flint::fmpz_set_form(tmp.d, (UWORD*)(symbol_stop+ABS(coeff_size/2)), ABS(coeff_size/2));
            // By construction, this is an exact division
            fmpz_divexact(coeff.d, coeff.d, tmp.d);
        }
 
        // Add the term to the poly
        fmpz_poly_set_coeff_fmpz(poly, exponent, coeff.d);
 
        if ( *term_stop == 0 ) {
            break;
        }
    }
 
    return arg_size;
}
/*
    #] flint::from_argument_poly :
 
    #[ flint::fmpz_get_form :
*/
// Write FORM's long integer representation of an fmpz at a, and put the number of WORDs at na.
// na carries the sign of the integer.
WORD flint::fmpz_get_form(fmpz_t z, WORD *a) {
 
    WORD na = 0;
    const int32_t sgn = fmpz_sgn(z);
    if ( sgn == -1 ) {
        fmpz_neg(z, z);
    }
    const int64_t nlimbs = fmpz_size(z);
 
    // This works but is UB?
    //fmpz_get_ui_array(reinterpret_cast<uint64_t*>(a), nlimbs, z);
 
    // Use fixed-size functions to get limb data where possible. These probably cover most real
    // cases. 
    if ( nlimbs == 1 ) {
        const uint64_t limb = fmpz_get_ui(z);
        a[0] = (WORD)(limb & 0xFFFFFFFF);
        na++;
        a[1] = (WORD)(limb >> BITSINWORD);
        if ( a[1] != 0 ) {
            na++;
        }
    }
    else if ( nlimbs == 2 ) {
        uint64_t limb_hi = 0, limb_lo = 0;
        fmpz_get_uiui((ulong*)&limb_hi, (ulong*)&limb_lo, z);
        a[0] = (WORD)(limb_lo & 0xFFFFFFFF);
            na++;
        a[1] = (WORD)(limb_lo >> BITSINWORD);
            na++;
        a[2] = (WORD)(limb_hi & 0xFFFFFFFF);
            na++;
        a[3] = (WORD)(limb_hi >> BITSINWORD);
        if ( a[3] != 0 ) {
            na++;
        }
    }
    else {
        vector<uint64_t> limb_data(nlimbs, 0);
        fmpz_get_ui_array((ulong*)limb_data.data(), (slong)nlimbs, z);
        for ( long i = 0; i < nlimbs; i++ ) {
            a[2*i] = (WORD)(limb_data[i] & 0xFFFFFFFF);
            na++;
            a[2*i+1] = (WORD)(limb_data[i] >> BITSINWORD);
            if ( a[2*i+1] != 0 || i < (nlimbs-1) ) {
                // The final limb might fit in a single 32bit WORD. Only
                // increment na if the final WORD is non zero.
                na++;
            }
        }
    }
 
    // And now put the sign in the number of limbs
    if ( sgn == -1 ) {
        na = -na;
    }
 
    return na;
}
/*
    #] flint::fmpz_get_form :
    #[ flint::fmpz_set_form :
*/
// Create an fmpz directly from FORM's long integer representation. fmpz uses 64bit unsigned limbs,
// but FORM uses 32bit UWORDs on 64bit architectures so we can't use fmpz_set_ui_array directly.
void flint::fmpz_set_form(fmpz_t z, UWORD *a, WORD na) {
 
    if ( na == 0 ) {
        fmpz_zero(z);
        return;
    }
 
    // Negative na represenents a negative number
    int32_t sgn = 1;
    if ( na < 0 ) {
        sgn = -1;
        na = -na;
    }
 
    // Remove padding. FORM stores numerators and denominators with equal numbers of limbs but we
    // don't need to do this within the fmpz. It is not necessary to do this really, the fmpz
    // doesn't add zero limbs unnecessarily, but we might be able to avoid creating the limb_data
    // array below.
    while ( a[na-1] == 0 ) {
        na--;
    }
 
    // If the number fits in fixed-size fmpz_set functions, we don't need to use additional memory
    // to convert to uint64_t. These probably cover most real cases.
    if ( na == 1 ) {
        fmpz_set_ui(z, (uint64_t)a[0]);
    }
    else if ( na == 2 ) {
        fmpz_set_ui(z, (((uint64_t)a[1])<<BITSINWORD) + (uint64_t)a[0]);
    }
    else if ( na == 3 ) {
        fmpz_set_uiui(z, (uint64_t)a[2], (((uint64_t)a[1])<<BITSINWORD) + (uint64_t)a[0]);
    }
    else if ( na == 4 ) {
        fmpz_set_uiui(z, (((uint64_t)a[3])<<BITSINWORD) + (uint64_t)a[2],
            (((uint64_t)a[1])<<BITSINWORD) + (uint64_t)a[0]);
    }
    else {
        const int32_t nlimbs = (na+1)/2;
        vector<uint64_t> limb_data(nlimbs, 0);
        for ( int32_t i = 0; i < nlimbs; i++ ) {
            if ( 2*i+1 <= na-1 ) {
                limb_data[i] = (uint64_t)a[2*i] + (((uint64_t)a[2*i+1])<<BITSINWORD);
            }
            else {
                limb_data[i] = (uint64_t)a[2*i];
            }
        }
        fmpz_set_ui_array(z, (ulong*)limb_data.data(), nlimbs);
    }
 
    // Finally set the sign.
    if ( sgn == -1 ) {
        fmpz_neg(z, z);
    }
 
    return;
}
/*
    #] flint::fmpz_set_form :
 
    #[ flint::gcd_mpoly :
*/
// Return a pointer to a buffer containing the GCD of the 0-terminated term lists at a and b.
// If must_fit_term, this should be a TermMalloc buffer. Otherwise Malloc1 the buffer.
// For multi-variate cases.
WORD* flint::gcd_mpoly(PHEAD const WORD *a, const WORD *b, const WORD must_fit_term,
    const var_map_t &var_map, const bool sort_vars) {
 
    flint::mpoly_ctx ctx(var_map.size());
    flint::mpoly pa(ctx.d), pb(ctx.d), denpa(ctx.d), denpb(ctx.d), gcd(ctx.d);
 
    flint::from_argument_mpoly(pa.d, denpa.d, a, false, var_map, ctx.d);
    flint::from_argument_mpoly(pb.d, denpb.d, b, false, var_map, ctx.d);
 
    // denpa, denpb must be 1:
    if ( fmpz_mpoly_is_one(denpa.d, ctx.d) != 1 ) {
        // INTERNAL_ERROR_EXCL_START
        MLOCK(ErrorMessageLock);
        MesPrint("!>flint::gcd_mpoly: error: denpa != 1");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
        // INTERNAL_ERROR_EXCL_STOP
    }
    if ( fmpz_mpoly_is_one(denpb.d, ctx.d) != 1 ) {
        // INTERNAL_ERROR_EXCL_START
        MLOCK(ErrorMessageLock);
        MesPrint("!>flint::gcd_mpoly: error: denpb != 1");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
        // INTERNAL_ERROR_EXCL_STOP
    }
 
    // poly returns pa if pa == pb, regardless of the lcoeff sign
    if ( fmpz_mpoly_equal(pa.d, pb.d, ctx.d) ) {
        fmpz_mpoly_set(gcd.d, pa.d, ctx.d);
    }
    else {
        // We need some gymnastics to have the same sign conventions as the poly class. It takes the
        // integer or univar content out of a,b, with the convention that the content sign matches
        // the lcoeff sign. Since FORM has already taken out the content, we are left with +-1. In
        // Flint, the content always has a positive sign so here we should find +1. Check this:
        flint::mpoly tmp(ctx.d);
        fmpz_mpoly_term_content(tmp.d, pa.d, ctx.d);
        if ( fmpz_mpoly_is_one(tmp.d, ctx.d) != 1 ) {
            // INTERNAL_ERROR_EXCL_START
            MLOCK(ErrorMessageLock);
            MesPrint("!>flint::gcd_mpoly: error: content of 1st arg != 1");
            MUNLOCK(ErrorMessageLock);
            Terminate(-1);
            // INTERNAL_ERROR_EXCL_STOP
        }
        fmpz_mpoly_term_content(tmp.d, pb.d, ctx.d);
        if ( fmpz_mpoly_is_one(tmp.d, ctx.d) != 1 ) {
            // INTERNAL_ERROR_EXCL_START
            MLOCK(ErrorMessageLock);
            MesPrint("!>flint::gcd_mpoly: error: content of 2nd arg != 1");
            MUNLOCK(ErrorMessageLock);
            Terminate(-1);
            // INTERNAL_ERROR_EXCL_STOP
        }
 
        // The poly class now divides the content out of a,b so that they have a positive lcoeff.
        // Then it multiplies the final gcd (which is given a positive lcoeff also) by
        // gcd(cont a, cont b). There it has gcd(1,1) = gcd(-1,1) = gcd(1,-1) = 1, and
        // gcd(-1,-1) = -1 (because of the pa==pb early return). So: if both input polys have a
        // negative lcoeff, we will flip the sign in the final result.
        bool flip_sign = 0;
        if ( ( fmpz_sgn(fmpz_mpoly_term_coeff_ref(pa.d, 0, ctx.d)) == -1 ) &&
            ( fmpz_sgn(fmpz_mpoly_term_coeff_ref(pb.d, 0, ctx.d)) == -1 ) ) {
            flip_sign = 1;
        }
 
        fmpz_mpoly_gcd(gcd.d, pa.d, pb.d, ctx.d);
        if ( flip_sign ) {
            fmpz_mpoly_neg(gcd.d, gcd.d, ctx.d);
        }
    }
 
    // This is freed by the caller
    WORD *res;
    if ( must_fit_term ) {
        res = TermMalloc("flint::gcd_mpoly");
    }
    else {
        // Determine the size of the GCD by passing write = false.
        const bool with_arghead = false;
        const bool write = false;
        const uint64_t prev_size = 0;
        const uint64_t gcd_size = (uint64_t)flint::to_argument_mpoly(BHEAD NULL,
            with_arghead, must_fit_term, write, prev_size, gcd.d, var_map, ctx.d, sort_vars);
        
        res = (WORD *)Malloc1(sizeof(WORD)*gcd_size, "flint::gcd_mpoly");
    }
 
    const bool with_arghead = false;
    const bool write = true;
    const uint64_t prev_size = 0;
    flint::to_argument_mpoly(BHEAD res, with_arghead, must_fit_term, write, prev_size, gcd.d,
        var_map, ctx.d, sort_vars);
 
    return res;
}
/*
    #] flint::gcd_mpoly :
    #[ flint::gcd_poly :
*/
// Return a pointer to a buffer containing the GCD of the 0-terminated term lists at a and b.
// If must_fit_term, this should be a TermMalloc buffer. Otherwise Malloc1 the buffer.
// For uni-variate cases.
WORD* flint::gcd_poly(PHEAD const WORD *a, const WORD *b, const WORD must_fit_term,
    const var_map_t &var_map) {
 
    flint::poly pa, pb, denpa, denpb, gcd;
 
    flint::from_argument_poly(pa.d, denpa.d, a, false);
    flint::from_argument_poly(pb.d, denpb.d, b, false);
 
    // denpa, denpb must be 1:
    if ( fmpz_poly_is_one(denpa.d) != 1 ) {
        // INTERNAL_ERROR_EXCL_START
        MLOCK(ErrorMessageLock);
        MesPrint("!>flint::gcd_poly: error: denpa != 1");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
        // INTERNAL_ERROR_EXCL_STOP
    }
    if ( fmpz_poly_is_one(denpb.d) != 1 ) {
        // INTERNAL_ERROR_EXCL_START
        MLOCK(ErrorMessageLock);
        MesPrint("!>flint::gcd_poly: error: denpb != 1");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
        // INTERNAL_ERROR_EXCL_STOP
    }
 
    // poly returns pa if pa == pb, regardless of the lcoeff sign
    if ( fmpz_poly_equal(pa.d, pb.d) ) {
        fmpz_poly_set(gcd.d, pa.d);
    }
    else {
        // Here, we don't have to make any sign flips like the mpoly case, because poly's
        // integer_gcd(1,1) = integer_gcd(-1,1) = integer_gcd(1,-1) = integer_gcd(-1,-1) = +1.
        // Still, verify that the content is 1:
        flint::fmpz tmp;
        fmpz_poly_content(tmp.d, pa.d);
        if ( fmpz_is_one(tmp.d) != 1 ) {
            // INTERNAL_ERROR_EXCL_START
            MLOCK(ErrorMessageLock);
            MesPrint("!>flint::gcd_poly: error: content of 1st arg != 1");
            MUNLOCK(ErrorMessageLock);
            Terminate(-1);
            // INTERNAL_ERROR_EXCL_STOP
        }
        fmpz_poly_content(tmp.d, pb.d);
        if ( fmpz_is_one(tmp.d) != 1 ) {
            // INTERNAL_ERROR_EXCL_START
            MLOCK(ErrorMessageLock);
            MesPrint("!>flint::gcd_poly: error: content of 2nd arg != 1");
            MUNLOCK(ErrorMessageLock);
            Terminate(-1);
            // INTERNAL_ERROR_EXCL_STOP
        }
 
        fmpz_poly_gcd(gcd.d, pa.d, pb.d);
    }
 
    // This is freed by the caller
    WORD *res;
    if ( must_fit_term ) {
        res = TermMalloc("flint::gcd_poly");
    }
    else {
        // Determine the size of the GCD by passing write = false.
        const bool with_arghead = false;
        const bool write = false;
        const uint64_t prev_size = 0;
        const uint64_t gcd_size = (uint64_t)flint::to_argument_poly(BHEAD NULL,
            with_arghead, must_fit_term, write, prev_size, gcd.d, var_map);
 
        res = (WORD *)Malloc1(sizeof(WORD)*gcd_size, "flint::gcd_poly");
    }
 
    const bool with_arghead = false;
    const uint64_t prev_size = 0;
    const bool write = true;
    flint::to_argument_poly(BHEAD res, with_arghead, must_fit_term, write, prev_size, gcd.d,
        var_map);
 
    return res;
}
/*
    #] flint::gcd_poly :
 
    #[ flint::get_variables :
*/
// Get the list of symbols which appear in the vector of expressions. These are polyratfun
// numerators, denominators or expressions from calls to gcd_ etc. Return this list as a map
// between indices and symbol codes.
// TODO FACTORSYMBOL last?
flint::var_map_t flint::get_variables(const vector <WORD *> &es, const bool with_arghead,
    const bool sort_vars) {
 
    int32_t num_vars = 0;
    // We count the total number of terms to determine "density".
    uint32_t num_terms = 0;
    // To be used if we sort by highest degree, as the poly code does.
    vector<int32_t> degrees;
    var_map_t var_map;
 
    // extract all variables
    for ( size_t ei = 0; ei < es.size(); ei++ ) {
        WORD *e = es[ei];
 
        // fast notation
        if ( *e == -SNUMBER ) {
            num_terms++;
        }
        else if ( *e == -SYMBOL ) {
            num_terms++;
            if ( !var_map.count(e[1]) ) {
                var_map[e[1]] = num_vars++;
                degrees.push_back(1);
            }
        }
        // JD: Here we need to check for non-symbol/number terms in fast notation.
        else if ( *e < 0 ) {
            MLOCK(ErrorMessageLock);
            MesPrint("ERROR: polynomials and polyratfuns must contain symbols only");
            MUNLOCK(ErrorMessageLock);
            Terminate(1);
        }
        else {
            for ( WORD i = with_arghead ? ARGHEAD:0; with_arghead ? i < e[0]:e[i] != 0; i += e[i] ) {
                num_terms++;
                const WORD coeff_size = e[i+e[i]-1];
                const WORD symbols_size = e[i] - ABS(coeff_size);
                if ( e[i+1] != SYMBOL && 1 < symbols_size ) {
                    MLOCK(ErrorMessageLock);
                    MesPrint("ERROR: polynomials and polyratfuns must contain symbols only");
                    MUNLOCK(ErrorMessageLock);
                    Terminate(1);
                }
 
                for ( WORD j = i+3; j < i+symbols_size; j += 2 ) {
                    const WORD symbol = e[j];
                    const WORD degree = e[j+1];
                    auto it = var_map.find(symbol);
                    if ( it == var_map.end() ) {
                        var_map[symbol] = num_vars++;
                        degrees.push_back(degree);
                    }
                    else {
                        degrees[it->second] = MaX(degrees[it->second], degree);
                    }
                }
            }
        }
    }
 
    if ( sort_vars ) {
        // bubble sort variables in decreasing order of degree
        // (this seems better for factorization)
        for ( size_t i = 0; i < var_map.size(); i++ ) {
            for ( size_t j = 0; j+1 < var_map.size(); j++ ) {
                if ( degrees[j] < degrees[j+1] ) {
                    swap(degrees[j], degrees[j+1]);
 
                    // Find the map keys associated with the values we want to swap
                    uint32_t j0 = 0;
                    uint32_t j1 = 0;
                    for ( auto x: var_map ) {
                        if ( x.second == j ) {
                            j0 = x.first;
                        }
                        else if ( x.second == j+1 ) {
                            j1 = x.first;
                        }
                    }
                    swap(var_map.at(j0), var_map.at(j1));
                }
            }
        }
    }
    // Otherwise, sort lexicographically in FORM's ordering
    else {
        for ( size_t i = 0; i < var_map.size(); i++ ) {
            for ( size_t j = 0; j+1 < var_map.size(); j++ ) {
                uint32_t j0 = 0;
                uint32_t j1 = 0;
                for ( auto x: var_map ) {
                    if ( x.second == j ) {
                        j0 = x.first;
                    }
                    else if ( x.second == j+1 ) {
                        j1 = x.first;
                    }
                }
                if ( j0 > j1 ) {
                    swap(var_map.at(j0), var_map.at(j1));
                }
            }
        }
    }
 
    if ( var_map.size() == 1 ) {
        // In the univariate case, if the polynomials are sufficiently sparse force the use of the
        // multivariate routines, which use a sparse representation, by adding a dummy map entry.
        if ( (float)num_terms <= UNIVARIATE_DENSITY_THR * (float)degrees[0] ) {
            // -1 will never be a symbol code. Built-in symbols from 0 to 19, and 20 is the first
            // user symbol.
            var_map[-1] = num_vars;
        }
    }
 
    return var_map;
}
/*
    #] flint::get_variables :
 
    #[ flint::inverse_poly :
*/
WORD* flint::inverse_poly(PHEAD const WORD *a, const WORD *b, const var_map_t &var_map) {
 
    flint::poly pa, pb, denpa, denpb;
 
    flint::from_argument_poly(pa.d, denpa.d, a, false);
    flint::from_argument_poly(pb.d, denpb.d, b, false);
 
    // fmpz_poly_xgcd is undefined if the content of pa and pb are not 1. Take the content out.
    // fmpz_poly_content gives the non-negative content and fmpz_poly_primitive normalizes to a
    // non-negative lcoeff, so we need to add the sign to the content if the polys have a negative
    // lcoeff. We don't need to keep the content of pb, it is a numerical multiple of the modulus.
    flint::fmpz content_a, resultant;
    flint::poly inverse, tmp;
    fmpz_poly_content(content_a.d, pa.d);
    if ( fmpz_sgn(fmpz_poly_lead(pa.d)) == -1 ) {
        fmpz_neg(content_a.d, content_a.d);
    }
    fmpz_poly_primitive_part(pa.d, pa.d);
    fmpz_poly_primitive_part(pb.d, pb.d);
 
    // Special cases:
    // Possibly strange that we give 1 for inverse_(x1,1) but here we take MMA's convention.
    if ( fmpz_poly_is_one(pa.d) && fmpz_poly_is_one(pb.d) ) {
        fmpz_poly_one(inverse.d);
        fmpz_one(resultant.d);
    }
    else if ( fmpz_poly_is_one(pb.d) ) {
        fmpz_poly_zero(inverse.d);
        fmpz_one(resultant.d);
    }
    else {
        // Now use the extended Euclidean algorithm to find inverse, resultant, tmp of the Bezout
        // identity: inverse*pa + tmp*pb = resultant. Then inverse/resultant is the multiplicative
        // inverse of pa mod pb. We'll divide by resultant in the denscale argument of to_argument_poly.
        fmpz_poly_xgcd(resultant.d, inverse.d, tmp.d, pa.d, pb.d);
 
        // If the resultant is zero, the inverse does not exist:
        if ( fmpz_is_zero(resultant.d) ) {
            MLOCK(ErrorMessageLock);
            MesPrint("flint::inverse_poly error: inverse does not exist");
            MUNLOCK(ErrorMessageLock);
            Terminate(-1);
        }
    }
 
    // Multiply inverse by denpa. denpb is a numerical multiple of the modulus, so doesn't matter.
    // We also need to divide by content_a, which we do in the denscale argument of to_argument_poly
    // by multiplying resultant by content_a here.
    fmpz_poly_mul(inverse.d, inverse.d, denpa.d);
    fmpz_mul(resultant.d, resultant.d, content_a.d);
 
 
    WORD* res;
    // First determine the result size, and malloc. The result should have no arghead. Here we use
    // the "scale" argument of to_argument_poly since resultant might not be 1.
    const bool with_arghead = false;
    bool write = false;
    const bool must_fit_term = false;
    const uint64_t prev_size = 0;
    const uint64_t res_size = (uint64_t)flint::to_argument_poly(BHEAD NULL,
        with_arghead, must_fit_term, write, prev_size, inverse.d, var_map, resultant.d);
    res = (WORD*)Malloc1(sizeof(WORD)*res_size, "flint::inverse_poly");
 
    write = true;
    flint::to_argument_poly(BHEAD res, with_arghead, must_fit_term, write, prev_size, inverse.d,
        var_map, resultant.d);
 
    return res;
}
/*
    #] flint::inverse_poly :
 
    #[ flint::mul_mpoly :
*/
// Return a pointer to a buffer containing the product of the 0-terminated term lists at a and b.
// For multi-variate cases.
WORD* flint::mul_mpoly(PHEAD const WORD *a, const WORD *b, const var_map_t &var_map,
    const bool sort_vars) {
 
    flint::mpoly_ctx ctx(var_map.size());
    flint::mpoly pa(ctx.d), pb(ctx.d), denpa(ctx.d), denpb(ctx.d);
 
    flint::from_argument_mpoly(pa.d, denpa.d, a, false, var_map, ctx.d);
    flint::from_argument_mpoly(pb.d, denpb.d, b, false, var_map, ctx.d);
 
    // denpa, denpb must be integers. Negative symbol powers have been converted to extra symbols.
    if ( fmpz_mpoly_is_fmpz(denpa.d, ctx.d) != 1 ) {
        // INTERNAL_ERROR_EXCL_START
        MLOCK(ErrorMessageLock);
        MesPrint("!>flint::mul_mpoly: error: denpa is non-constant");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
        // INTERNAL_ERROR_EXCL_STOP
    }
    if ( fmpz_mpoly_is_fmpz(denpb.d, ctx.d) != 1 ) {
        // INTERNAL_ERROR_EXCL_START
        MLOCK(ErrorMessageLock);
        MesPrint("!>flint::mul_mpoly: error: denpb is non-constant");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
        // INTERNAL_ERROR_EXCL_STOP
    }
 
    // Multiply numerators, store result in pa
    fmpz_mpoly_mul(pa.d, pa.d, pb.d, ctx.d);
    // Multiply denominators, store result in denpa, and convert to an fmpz:
    fmpz_mpoly_mul(denpa.d, denpa.d, denpb.d, ctx.d);
    flint::fmpz den;
    fmpz_mpoly_get_fmpz(den.d, denpa.d, ctx.d);
 
    WORD* res;
    // First determine the result size, and malloc. The result should have no arghead. Here we use
    // the "scale" argument of to_argument_mpoly since den might not be 1.
    const bool with_arghead = false;
    bool write = false;
    const bool must_fit_term = false;
    const uint64_t prev_size = 0;
    const uint64_t mul_size = (uint64_t)flint::to_argument_mpoly(BHEAD NULL,
        with_arghead, must_fit_term, write, prev_size, pa.d, var_map, ctx.d, sort_vars, den.d);
    res = (WORD*)Malloc1(sizeof(WORD)*mul_size, "flint::mul_mpoly");
 
    write = true;
    flint::to_argument_mpoly(BHEAD res, with_arghead, must_fit_term, write, prev_size, pa.d,
        var_map, ctx.d, sort_vars, den.d);
 
    return res;
}
/*
    #] flint::mul_mpoly :
    #[ flint::mul_poly :
*/
// Return a pointer to a buffer containing the product of the 0-terminated term lists at a and b.
// For uni-variate cases.
WORD* flint::mul_poly(PHEAD const WORD *a, const WORD *b, const var_map_t &var_map) {
 
    flint::poly pa, pb, denpa, denpb;
 
    flint::from_argument_poly(pa.d, denpa.d, a, false);
    flint::from_argument_poly(pb.d, denpb.d, b, false);
 
    // denpa, denpb must be integers. Negative symbol powers have been converted to extra symbols.
    if ( fmpz_poly_degree(denpa.d) != 0 ) {
        // INTERNAL_ERROR_EXCL_START
        MLOCK(ErrorMessageLock);
        MesPrint("!>flint::mul_poly: error: denpa is non-constant");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
        // INTERNAL_ERROR_EXCL_STOP
    }
    if ( fmpz_poly_degree(denpb.d) != 0 ) {
        // INTERNAL_ERROR_EXCL_START
        MLOCK(ErrorMessageLock);
        MesPrint("!>flint::mul_poly: error: denpb is non-constant");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
        // INTERNAL_ERROR_EXCL_STOP
    }
 
    // Multiply numerators, store result in pa
    fmpz_poly_mul(pa.d, pa.d, pb.d);
    // Multiply denominators, store result in denpa, and convert to an fmpz:
    fmpz_poly_mul(denpa.d, denpa.d, denpb.d);
    flint::fmpz den;
    fmpz_poly_get_coeff_fmpz(den.d, denpa.d, 0);
 
    WORD* res;
    // First determine the result size, and malloc. The result should have no arghead. Here we use
    // the "scale" argument of to_argument_poly since den might not be 1.
    const bool with_arghead = false;
    bool write = false;
    const bool must_fit_term = false;
    const uint64_t prev_size = 0;
    const uint64_t mul_size = (uint64_t)flint::to_argument_poly(BHEAD NULL,
        with_arghead, must_fit_term, write, prev_size, pa.d, var_map, den.d);
    res = (WORD*)Malloc1(sizeof(WORD)*mul_size, "flint::mul_poly");
 
    write = true;
    flint::to_argument_poly(BHEAD res, with_arghead, must_fit_term, write, prev_size, pa.d,
        var_map, den.d);
 
    return res;
}
/*
    #] flint::mul_poly :
 
    #[ flint::ratfun_add_mpoly :
*/
// Add the multi-variate FORM rational polynomials at t1 and t2. The result is written at out.
void flint::ratfun_add_mpoly(PHEAD const WORD *t1, const WORD *t2, WORD *out,
    const var_map_t &var_map, const bool sort_vars) {
 
    flint::mpoly_ctx ctx(var_map.size());
    flint::mpoly gcd(ctx.d), num1(ctx.d), den1(ctx.d), num2(ctx.d), den2(ctx.d);
 
    flint::ratfun_read_mpoly(t1, num1.d, den1.d, var_map, ctx.d);
    flint::ratfun_read_mpoly(t2, num2.d, den2.d, var_map, ctx.d);
 
    if ( fmpz_mpoly_cmp(den1.d, den2.d, ctx.d) != 0 ) {
        fmpz_mpoly_gcd_cofactors(gcd.d, den1.d, den2.d, den1.d, den2.d, ctx.d);
 
        fmpz_mpoly_mul(num1.d, num1.d, den2.d, ctx.d);
        fmpz_mpoly_mul(num2.d, num2.d, den1.d, ctx.d);
 
        fmpz_mpoly_add(num1.d, num1.d, num2.d, ctx.d);
        fmpz_mpoly_mul(den1.d, den1.d, den2.d, ctx.d);
        fmpz_mpoly_mul(den1.d, den1.d, gcd.d,  ctx.d);
    }
    else {
        fmpz_mpoly_add(num1.d, num1.d, num2.d, ctx.d);
    }
 
    // Finally divide out any common factors between the resulting num1, den1:
    fmpz_mpoly_gcd_cofactors(gcd.d, num1.d, den1.d, num1.d, den1.d, ctx.d);
 
    flint::util::fix_sign_fmpz_mpoly_ratfun(num1.d, den1.d, ctx.d);
 
    // Result in FORM notation:
    *out++ = AR.PolyFun;
    WORD* args_size = out++;
    WORD* args_flag = out++;
    *args_flag = 0; // clean prf
    FILLFUN3(out); // Remainder of funhead, if it is larger than 3
 
    const bool with_arghead = true;
    const bool must_fit_term = true;
    const bool write = true;
    // prev_size + 4, to account for final term size and coeff of "1/1"
    out += flint::to_argument_mpoly(BHEAD out, with_arghead, must_fit_term, write, out-args_size+4,
        num1.d, var_map, ctx.d, sort_vars);
    out += flint::to_argument_mpoly(BHEAD out, with_arghead, must_fit_term, write, out-args_size+4,
        den1.d, var_map, ctx.d, sort_vars);
 
    *args_size = out - args_size + 1; // The +1 is to include the function ID
    AT.WorkPointer = out;
}
/*
    #] flint::ratfun_add_mpoly :
    #[ flint::ratfun_add_poly :
*/
// Add the uni-variate FORM rational polynomials at t1 and t2. The result is written at out.
void flint::ratfun_add_poly(PHEAD const WORD *t1, const WORD *t2, WORD *out,
    const var_map_t &var_map) {
 
    flint::poly gcd, num1, den1, num2, den2;
 
    flint::ratfun_read_poly(t1, num1.d, den1.d);
    flint::ratfun_read_poly(t2, num2.d, den2.d);
 
    if ( fmpz_poly_equal(den1.d, den2.d) == 0 ) {
        flint::util::simplify_fmpz_poly(den1.d, den2.d, gcd.d);
 
        fmpz_poly_mul(num1.d, num1.d, den2.d);
        fmpz_poly_mul(num2.d, num2.d, den1.d);
 
        fmpz_poly_add(num1.d, num1.d, num2.d);
        fmpz_poly_mul(den1.d, den1.d, den2.d);
        fmpz_poly_mul(den1.d, den1.d, gcd.d);
    }
    else {
        fmpz_poly_add(num1.d, num1.d, num2.d);
    }
 
    // Finally divide out any common factors between the resulting num1, den1:
    flint::util::simplify_fmpz_poly(num1.d, den1.d, gcd.d);
 
    flint::util::fix_sign_fmpz_poly_ratfun(num1.d, den1.d);
 
    // Result in FORM notation:
    *out++ = AR.PolyFun;
    WORD* args_size = out++;
    WORD* args_flag = out++;
    *args_flag = 0; // clean prf
    FILLFUN3(out); // Remainder of funhead, if it is larger than 3
 
    const bool with_arghead = true;
    const bool must_fit_term = true;
    const bool write = true;
    // prev_size + 4, to account for final term size and coeff of "1/1"
    out += flint::to_argument_poly(BHEAD out, with_arghead, must_fit_term, write, out-args_size+4,
        num1.d, var_map);
    out += flint::to_argument_poly(BHEAD out, with_arghead, must_fit_term, write, out-args_size+4,
        den1.d, var_map);
 
    *args_size = out - args_size + 1; // The +1 is to include the function ID
    AT.WorkPointer = out;
}
/*
    #] flint::ratfun_add_poly :
 
    #[ flint::ratfun_normalize_mpoly :
*/
// Multiply and simplify occurrences of the multi-variate FORM rational polynomials found in term.
// The final term is written in place, with the rational polynomial at the end.
void flint::ratfun_normalize_mpoly(PHEAD WORD *term, const var_map_t &var_map,
    const bool sort_vars) {
 
    // The length of the coefficient
    const WORD ncoeff = (term + *term)[-1];
    // The end of the term data, before the coefficient:
    const WORD *tstop = term + *term - ABS(ncoeff);
 
    flint::mpoly_ctx ctx(var_map.size());
    flint::mpoly num1(ctx.d), den1(ctx.d), num2(ctx.d), den2(ctx.d), gcd(ctx.d);
 
    // Start with "trivial" polynomials, and multiply in the num and den of the prf which appear.
    flint::fmpz tmpNum, tmpDen;
    flint::fmpz_set_form(tmpNum.d, (UWORD*)tstop, ncoeff/2);
    flint::fmpz_set_form(tmpDen.d, (UWORD*)tstop+ABS(ncoeff/2), ABS(ncoeff/2));
    fmpz_mpoly_set_fmpz(num1.d, tmpNum.d, ctx.d);
    fmpz_mpoly_set_fmpz(den1.d, tmpDen.d, ctx.d);
 
    // Loop over the occurrences of PolyFun in the term, and multiply in to num1, den1.
    // s tracks where we are writing the non-PolyFun term data. The final PolyFun will
    // go at the end.
    WORD* term_size = term;
    WORD* s = term + 1;
    for ( WORD *t = term + 1; t < tstop; ) {
        if ( *t == AR.PolyFun ) {
            flint::ratfun_read_mpoly(t, num2.d, den2.d, var_map, ctx.d);
 
            // get gcd of num1,den2 and num2,den1 and then assemble
            fmpz_mpoly_gcd_cofactors(gcd.d, num1.d, den2.d, num1.d, den2.d, ctx.d);
            fmpz_mpoly_gcd_cofactors(gcd.d, num2.d, den1.d, num2.d, den1.d, ctx.d);
 
            fmpz_mpoly_mul(num1.d, num1.d, num2.d, ctx.d);
            fmpz_mpoly_mul(den1.d, den1.d, den2.d, ctx.d);
 
            t += t[1];
        }
 
        else {
            // Not a PolyFun, just copy or skip over
            WORD i = t[1];
            if ( s != t ) { NCOPY(s,t,i); }
            else { t += i; s += i; }
        }
    }
 
    flint::util::fix_sign_fmpz_mpoly_ratfun(num1.d, den1.d, ctx.d);
 
    // Result in FORM notation:
    WORD* out = s;
    *out++ = AR.PolyFun;
    WORD* args_size = out++;
    WORD* args_flag = out++;
    *args_flag &= ~MUSTCLEANPRF;
 
    const bool with_arghead = true;
    const bool must_fit_term = true;
    const bool write = true;
    out += flint::to_argument_mpoly(BHEAD out, with_arghead, must_fit_term, write, out-term_size,
        num1.d, var_map, ctx.d, sort_vars);
    out += flint::to_argument_mpoly(BHEAD out, with_arghead, must_fit_term, write, out-term_size,
        den1.d, var_map, ctx.d, sort_vars);
 
    *args_size = out - args_size + 1; // The +1 is to include the function ID
 
    // +3 for the coefficient of 1/1, which is added after the check
    if ( sizeof(WORD)*(out-term_size+3) > (size_t)AM.MaxTer ) {
        MLOCK(ErrorMessageLock);
        MesPrint("flint::ratfun_normalize: output exceeds MaxTermSize");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
    }
 
    *out++ = 1;
    *out++ = 1;
    *out++ = 3; // the term's coefficient is now 1/1
 
    *term_size = out - term_size;
}
/*
    #] flint::ratfun_normalize_mpoly :
    #[ flint::ratfun_normalize_poly :
*/
// Multiply and simplify occurrences of the uni-variate FORM rational polynomials found in term.
// The final term is written in place, with the rational polynomial at the end.
void flint::ratfun_normalize_poly(PHEAD WORD *term, const var_map_t &var_map) {
 
    // The length of the coefficient
    const WORD ncoeff = (term + *term)[-1];
    // The end of the term data, before the coefficient:
    const WORD *tstop = term + *term - ABS(ncoeff);
 
    flint::poly num1, den1, num2, den2, gcd;
 
    // Start with "trivial" polynomials, and multiply in the num and den of the prf which appear.
    flint::fmpz tmpNum, tmpDen;
    flint::fmpz_set_form(tmpNum.d, (UWORD*)tstop, ncoeff/2);
    flint::fmpz_set_form(tmpDen.d, (UWORD*)tstop+ABS(ncoeff/2), ABS(ncoeff/2));
    fmpz_poly_set_fmpz(num1.d, tmpNum.d);
    fmpz_poly_set_fmpz(den1.d, tmpDen.d);
 
    // Loop over the occurrences of PolyFun in the term, and multiply in to num1, den1.
    // s tracks where we are writing the non-PolyFun term data. The final PolyFun will
    // go at the end.
    WORD* term_size = term;
    WORD* s = term + 1;
    for ( WORD *t = term + 1; t < tstop; ) {
        if ( *t == AR.PolyFun ) {
            flint::ratfun_read_poly(t, num2.d, den2.d);
 
            // get gcd of num1,den2 and num2,den1 and then assemble
            flint::util::simplify_fmpz_poly(num1.d, den2.d, gcd.d);
            flint::util::simplify_fmpz_poly(num2.d, den1.d, gcd.d);
 
            fmpz_poly_mul(num1.d, num1.d, num2.d);
            fmpz_poly_mul(den1.d, den1.d, den2.d);
 
            t += t[1];
        }
 
        else {
            // Not a PolyFun, just copy or skip over
            WORD i = t[1];
            if ( s != t ) { NCOPY(s,t,i); }
            else { t += i; s += i; }
        }
    }
 
    flint::util::fix_sign_fmpz_poly_ratfun(num1.d, den1.d);
 
    // Result in FORM notation:
    WORD* out = s;
    *out++ = AR.PolyFun;
    WORD* args_size = out++;
    WORD* args_flag = out++;
    *args_flag &= ~MUSTCLEANPRF;
 
    const bool with_arghead = true;
    const bool must_fit_term = true;
    const bool write = true;
    out += flint::to_argument_poly(BHEAD out, with_arghead, must_fit_term, write, out-term_size,
        num1.d, var_map);
    out += flint::to_argument_poly(BHEAD out, with_arghead, must_fit_term, write, out-term_size,
        den1.d, var_map);
 
    *args_size = out - args_size + 1; // The +1 is to include the function ID
 
    // +3 for the coefficient of 1/1, which is added after the check
    if ( sizeof(WORD)*(out-term_size+3) > (size_t)AM.MaxTer ) {
        MLOCK(ErrorMessageLock);
        MesPrint("flint::ratfun_normalize: output exceeds MaxTermSize");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
    }
 
    *out++ = 1;
    *out++ = 1;
    *out++ = 3; // the term's coefficient is now 1/1
 
    *term_size = out - term_size;
}
/*
    #] flint::ratfun_normalize_poly :
 
    #[ flint::ratfun_read_mpoly :
*/
// Read the multi-variate FORM rational polynomial at a and create fmpz_mpoly_t numerator and
// denominator.
void flint::ratfun_read_mpoly(const WORD *a, fmpz_mpoly_t num, fmpz_mpoly_t den,
    const var_map_t &var_map, fmpz_mpoly_ctx_t ctx) {
 
    // The end of the arguments:
    const WORD* arg_stop = a+a[1];
 
    const bool must_normalize = (a[2] & MUSTCLEANPRF) != 0;
 
    a += FUNHEAD;
    if ( a >= arg_stop ) {
        MLOCK(ErrorMessageLock);
        MesPrint("ERROR: PolyRatFun cannot have zero arguments");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
    }
 
    // Polys to collect the "den of the num" and "den of the den".
    // Input can arrive like this when enabling the PolyRatFun or moving things into it.
    flint::mpoly den_num(ctx), den_den(ctx);
 
    // Read the numerator
    flint::from_argument_mpoly(num, den_num.d, a, true, var_map, ctx);
    NEXTARG(a);
 
    if ( a < arg_stop ) {
        // Read the denominator
        flint::from_argument_mpoly(den, den_den.d, a, true, var_map, ctx);
        NEXTARG(a);
    }
    else {
        // The denominator is 1
/* INTERNAL_ERROR_EXCL_START */
        MLOCK(ErrorMessageLock);
        MesPrint("!>implement this");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
/* INTERNAL_ERROR_EXCL_STOP */
    }
    if ( a < arg_stop ) {
        MLOCK(ErrorMessageLock);
        MesPrint("ERROR: PolyRatFun cannot have more than two arguments");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
    }
 
    // Multiply the num by den_den and den by den_num:
    fmpz_mpoly_mul(num, num, den_den.d, ctx);
    fmpz_mpoly_mul(den, den, den_num.d, ctx);
 
    if ( must_normalize ) {
        flint::mpoly gcd(ctx);
        fmpz_mpoly_gcd_cofactors(gcd.d, num, den, num, den, ctx);
    }
}
/*
    #] flint::ratfun_read_mpoly :
    #[ flint::ratfun_read_poly :
*/
// Read the uni-variate FORM rational polynomial at a and create fmpz_mpoly_t numerator and
// denominator.
void flint::ratfun_read_poly(const WORD *a, fmpz_poly_t num, fmpz_poly_t den) {
 
    // The end of the arguments:
    const WORD* arg_stop = a+a[1];
 
    const bool must_normalize = (a[2] & MUSTCLEANPRF) != 0;
 
    a += FUNHEAD;
    if ( a >= arg_stop ) {
        MLOCK(ErrorMessageLock);
        MesPrint("ERROR: PolyRatFun cannot have zero arguments");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
    }
 
    // Polys to collect the "den of the num" and "den of the den".
    // Input can arrive like this when enabling the PolyRatFun or moving things into it.
    flint::poly den_num, den_den;
 
    // Read the numerator
    flint::from_argument_poly(num, den_num.d, a, true);
    NEXTARG(a);
 
    if ( a < arg_stop ) {
        // Read the denominator
        flint::from_argument_poly(den, den_den.d, a, true);
        NEXTARG(a);
    }
    else {
        // The denominator is 1
/* INTERNAL_ERROR_EXCL_START */
        MLOCK(ErrorMessageLock);
        MesPrint("!>implement this");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
/* INTERNAL_ERROR_EXCL_STOP */
    }
    if ( a < arg_stop ) {
        MLOCK(ErrorMessageLock);
        MesPrint("ERROR: PolyRatFun cannot have more than two arguments");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
    }
 
    // Multiply the num by den_den and den by den_num:
    fmpz_poly_mul(num, num, den_den.d);
    fmpz_poly_mul(den, den, den_num.d);
 
    if ( must_normalize ) {
        flint::poly gcd;
        flint::util::simplify_fmpz_poly(num, den, gcd.d);
    }
}
/*
    #] flint::ratfun_read_poly :
    #[ flint::to_argument_mpoly :
*/
// Convert a fmpz_mpoly_t to a FORM argument (or 0-terminated list of terms: with_arghead==false).
// If the caller is building an output term, prev_size contains the size of the term so far, to
// check that the output fits in AM.MaxTer if must_fit_term.
// All coefficients will be divided by denscale (which might just be 1).
// If write is false, we never write to out but only track the total would-be size. This lets this
// function be repurposed as a "size of FORM notation" function without duplicating the code.
#define IFW(x) { if ( write ) {x;} }
uint64_t flint::to_argument_mpoly(PHEAD WORD *out, const bool with_arghead,
    const bool must_fit_term, const bool write, const uint64_t prev_size, const fmpz_mpoly_t poly,
    const var_map_t &var_map, const fmpz_mpoly_ctx_t ctx, const bool sort_vars,
    const fmpz_t denscale) {
 
    // out is modified later, keep the pointer at entry
    const WORD* out_entry = out;
 
    // Track the total size written. We could do this with pointer differences, but if
    // write == false we don't write to or move out to be able to find the size that way.
    uint64_t ws = 0;
 
    // Check there is at least space for ARGHEAD WORDs (the arghead or fast-notation number/symbol)
    if ( write && must_fit_term && (sizeof(WORD)*(prev_size + ARGHEAD) > (size_t)AM.MaxTer) ) {
        MLOCK(ErrorMessageLock);
        MesPrint("flint::to_argument_mpoly: output exceeds MaxTermSize");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
    }
 
    // Create the inverse of var_map, so we don't have to search it for each symbol written
    vector<uint32_t> var_map_inv;
    var_map_inv.resize(var_map.size());
    for ( auto x : var_map ) {
        var_map_inv[x.second] = x.first;
    }
 
    vector<int64_t> exponents(var_map.size());
    const int64_t n_terms = fmpz_mpoly_length(poly, ctx);
 
    if ( n_terms == 0 ) {
        if ( with_arghead ) {
            IFW(*out++ = -SNUMBER); ws++;
            IFW(*out++ = 0); ws++;
            return ws;
        }
        else {
            IFW(*out++ = 0); ws++;
            return ws;
        }
    }
 
    // For dividing out denscale
    flint::fmpz coeff, den, gcd;
 
    // The mpoly might be constant or a single symbol with coeff 1. Use fast notation if possible.
    if ( with_arghead && n_terms == 1 ) {
 
        if ( fmpz_mpoly_is_fmpz(poly, ctx) ) {
            // The mpoly is constant. Use fast notation if the number is an integer and small enough:
 
            fmpz_mpoly_get_term_coeff_fmpz(coeff.d, poly, 0, ctx);
            fmpz_set(den.d, denscale);
            flint::util::simplify_fmpz(coeff.d, den.d, gcd.d);
 
            if ( fmpz_is_one(den.d) && fmpz_fits_si(coeff.d) ) {
                const int64_t fast_coeff = fmpz_get_si(coeff.d);
                // While ">=", could work here, FORM does not use fast notation for INT_MIN
                if ( fast_coeff > INT32_MIN && fast_coeff <= INT32_MAX ) {
                    IFW(*out++ = -SNUMBER); ws++;
                    IFW(*out++ = (WORD)fast_coeff); ws++;
                    return ws;
                }
            }
        }
 
        else {
            fmpz_mpoly_get_term_coeff_fmpz(coeff.d, poly, 0, ctx);
            fmpz_set(den.d, denscale);
            flint::util::simplify_fmpz(coeff.d, den.d, gcd.d);
 
            if ( fmpz_is_one(coeff.d) && fmpz_is_one(den.d) ) {
                // The coefficient is one. Now check the symbol powers:
 
                fmpz_mpoly_get_term_exp_si((slong*)exponents.data(), poly, 0, ctx);
                int64_t use_fast = 0;
                uint32_t fast_symbol = 0;
 
                for ( size_t i = 0; i < var_map.size(); i++ ) {
                    if ( exponents[i] == 1 ) fast_symbol = var_map_inv[i];
                    use_fast += exponents[i];
                }
 
                // use_fast has collected the total degree. If it is 1, then fast_symbol holds the code
                if ( use_fast == 1 ) {
                    IFW(*out++ = -SYMBOL); ws++;
                    IFW(*out++ = fast_symbol); ws++;
                    return ws;
                }
            }
        }
    }
 
 
    WORD *tmp_coeff = (WORD *)NumberMalloc("flint::to_argument_mpoly");
    WORD *tmp_den = (WORD *)NumberMalloc("flint::to_argument_mpoly");
 
    WORD* arg_size = 0;
    WORD* arg_flag = 0;
    if ( with_arghead ) {
        IFW(arg_size = out++); ws++; // total arg size
        IFW(arg_flag = out++); ws++;
        IFW(*arg_flag = 0); // clean argument
    }
 
    for ( int64_t i = 0; i < n_terms; i++ ) {
 
        fmpz_mpoly_get_term_exp_si((slong*)exponents.data(), poly, i, ctx);
 
        fmpz_mpoly_get_term_coeff_fmpz(coeff.d, poly, i, ctx);
        fmpz_set(den.d, denscale);
        flint::util::simplify_fmpz(coeff.d, den.d, gcd.d);
 
        uint32_t num_symbols = 0;
        for ( size_t j = 0; j < var_map.size(); j++ ) {
            if ( exponents[j] != 0 ) { num_symbols += 1; }
        }
 
        // Convert the coefficient, write in temporary space
        const WORD num_size = flint::fmpz_get_form(coeff.d, tmp_coeff);
        const WORD den_size = flint::fmpz_get_form(den.d, tmp_den);
        const WORD coeff_size = [num_size, den_size] () -> WORD {
            WORD size = ABS(num_size) > ABS(den_size) ? ABS(num_size) : ABS(den_size);
            return size * SGN(num_size) * SGN(den_size);
        }();
 
        // Now we have the number of symbols and the coeff size, we can determine the output size.
        // Check it fits if necessary: term size, num,den of "coeff_size", +1 for total coeff size
        uint64_t current_size = prev_size + ws + 1 + 2*ABS(coeff_size) + 1;
        if ( num_symbols ) {
            // symbols header, code,power of each symbol:
            current_size += 2 + 2*num_symbols;
        }
        if ( write && must_fit_term && (sizeof(WORD)*current_size > (size_t)AM.MaxTer) ) {
            MLOCK(ErrorMessageLock);
            MesPrint("flint::to_argument_mpoly: output exceeds MaxTermSize");
            MUNLOCK(ErrorMessageLock);
            Terminate(-1);
        }
 
        WORD* term_size = 0;
        IFW(term_size = out++); ws++;
        if ( num_symbols ) {
            IFW(*out++ = SYMBOL); ws++;
            WORD* symbol_size = 0;
            IFW(symbol_size = out++); ws++;
            IFW(*symbol_size = 2);
 
            for ( size_t j = 0; j < var_map.size(); j++ ) {
                if ( exponents[j] != 0 ) {
                    IFW(*out++ = var_map_inv[j]); ws++;
                    IFW(*out++ = exponents[j]); ws++;
                    IFW(*symbol_size += 2);
                }
            }
        }
 
        // Copy numerator
        for ( WORD j = 0; j < ABS(num_size); j++ ) {
            IFW(*out++ = tmp_coeff[j]); ws++;
        }
        for ( WORD j = ABS(num_size); j < ABS(coeff_size); j++ ) {
            IFW(*out++ = 0); ws++;
        }
        // Copy denominator
        for ( WORD j = 0; j < ABS(den_size); j++ ) {
            IFW(*out++ = tmp_den[j]); ws++;
        }
        for ( WORD j = ABS(den_size); j < ABS(coeff_size); j++ ) {
            IFW(*out++ = 0); ws++;
        }
 
        IFW(*out = 2*ABS(coeff_size) + 1); // the size of the coefficient
        IFW(if ( coeff_size < 0 ) { *out = -(*out); });
        IFW(out++); ws++;
 
        IFW(*term_size = out - term_size);
    }
 
    if ( with_arghead ) {
        IFW(*arg_size = out - arg_size);
        if ( write && sort_vars ) {
            // Sort into form highfirst ordering, if we have potentially re-ordered the variables by
            // degree, in flint::get_variables. Otherwise, we can rely on FLINT's sorting with
            // ORD_LEX and the variables ordered in FORM's lexicographic order.
            flint::form_sort(BHEAD (WORD*)(out_entry));
        }
    }
    else {
        // with no arghead, we write a terminating zero
        IFW(*out++ = 0); ws++;
    }
 
    NumberFree(tmp_coeff, "flint::to_argument_mpoly");
    NumberFree(tmp_den, "flint::to_argument_mpoly");
 
    return ws;
}
 
// If no denscale argument is supplied, just set it to 1 and call the usual function
uint64_t flint::to_argument_mpoly(PHEAD WORD *out, const bool with_arghead,
    const bool must_fit_term, const bool write, const uint64_t prev_size, const fmpz_mpoly_t poly,
    const var_map_t &var_map, const fmpz_mpoly_ctx_t ctx, const bool sort_vars) {
 
    flint::fmpz tmp;
    fmpz_set_ui(tmp.d, 1);
 
    uint64_t ret = flint::to_argument_mpoly(BHEAD out, with_arghead, must_fit_term, write,
        prev_size, poly, var_map, ctx, sort_vars, tmp.d);
 
    return ret;
}
/*
    #] flint::to_argument_mpoly :
    #[ flint::to_argument_poly :
*/
// Convert a fmpz_poly_t to a FORM argument (or 0-terminated list of terms: with_arghead==false).
// If the caller is building an output term, prev_size contains the size of the term so far, to
// check that the output fits in AM.MaxTer if must_fit_term.
// All coefficients will be divided by denscale (which might just be 1).
// If write is false, we never write to out but only track the total would-be size. This lets this
// function be repurposed as a "size of FORM notation" function without duplicating the code.
uint64_t flint::to_argument_poly(PHEAD WORD *out, const bool with_arghead,
    const bool must_fit_term, const bool write, const uint64_t prev_size, const fmpz_poly_t poly,
    const var_map_t &var_map, const fmpz_t denscale) {
 
    // Track the total size written. We could do this with pointer differences, but if
    // write == false we don't write to or move out to be able to find the size that way.
    uint64_t ws = 0;
 
    // Check there is at least space for ARGHEAD WORDs (the arghead or fast-notation number/symbol)
    if ( write && must_fit_term && (sizeof(WORD)*(prev_size + ARGHEAD) > (size_t)AM.MaxTer) ) {
        MLOCK(ErrorMessageLock);
        MesPrint("flint::to_argument_poly: output exceeds MaxTermSize");
        MUNLOCK(ErrorMessageLock);
        Terminate(-1);
    }
 
    // Create the inverse of var_map, so we don't have to search it for each symbol written
    vector<uint32_t> var_map_inv;
    var_map_inv.resize(var_map.size());
    for ( auto x : var_map ) {
        var_map_inv[x.second] = x.first;
    }
 
    const int64_t n_terms = fmpz_poly_length(poly);
 
    // The poly is zero
    if ( n_terms == 0 ) {
        if ( with_arghead ) {
            IFW(*out++ = -SNUMBER); ws++;
            IFW(*out++ = 0); ws++;
            return ws;
        }
        else {
            IFW(*out++ = 0); ws++;
            return ws;
        }
    }
 
    // For dividing out denscale
    flint::fmpz coeff, den, gcd;
 
    // The poly is constant, use fast notation if the coefficient is integer and small enough
    if ( with_arghead && n_terms == 1 ) {
 
        fmpz_poly_get_coeff_fmpz(coeff.d, poly, 0);
        fmpz_set(den.d, denscale);
        flint::util::simplify_fmpz(coeff.d, den.d, gcd.d);
 
        if ( fmpz_is_one(den.d) && fmpz_fits_si(coeff.d) ) {
            const long fast_coeff = fmpz_get_si(coeff.d);
            // While ">=", could work here, FORM does not use fast notation for INT_MIN
            if ( fast_coeff > INT_MIN && fast_coeff <= INT_MAX ) {
                IFW(*out++ = -SNUMBER); ws++;
                IFW(*out++ = (WORD)fast_coeff); ws++;
                return ws;
            }
        }
    }
 
    // The poly might be a single symbol with coeff 1, use fast notation if so.
    if ( with_arghead && n_terms == 2 ) {
        if ( fmpz_is_zero(fmpz_poly_get_coeff_ptr(poly, 0)) ) {
            // The constant term is zero
 
            fmpz_poly_get_coeff_fmpz(coeff.d, poly, 1);
            fmpz_set(den.d, denscale);
            flint::util::simplify_fmpz(coeff.d, den.d, gcd.d);
 
            if ( fmpz_is_one(coeff.d) && fmpz_is_one(den.d) ) {
                // Single symbol with coeff 1. Use fast notation:
                IFW(*out++ = -SYMBOL); ws++;
                IFW(*out++ = var_map_inv[0]); ws++;
                return ws;
            }
        }
    }
 
    WORD *tmp_coeff = (WORD *)NumberMalloc("flint::to_argument_poly");
    WORD *tmp_den = (WORD *)NumberMalloc("flint::to_argument_mpoly");
 
    WORD* arg_size = 0;
    WORD* arg_flag = 0;
    if ( with_arghead ) {
        IFW(arg_size = out++); ws++; // total arg size
        IFW(arg_flag = out++); ws++;
        IFW(*arg_flag = 0); // clean argument
    }
 
    // In reverse, since we want a "highfirst" output
    for ( int64_t i = n_terms-1; i >= 0; i-- ) {
 
        // fmpz_poly is dense, there might be many zero coefficients:
        if ( !fmpz_is_zero(fmpz_poly_get_coeff_ptr(poly, i)) ) {
 
            fmpz_poly_get_coeff_fmpz(coeff.d, poly, i);
            fmpz_set(den.d, denscale);
            flint::util::simplify_fmpz(coeff.d, den.d, gcd.d);
 
            // Convert the coefficient, write in temporary space
            const WORD num_size = flint::fmpz_get_form(coeff.d, tmp_coeff);
            const WORD den_size = flint::fmpz_get_form(den.d, tmp_den);
            const WORD coeff_size = [num_size, den_size] () -> WORD {
                WORD size = ABS(num_size) > ABS(den_size) ? ABS(num_size) : ABS(den_size);
                return size * SGN(num_size) * SGN(den_size);
            }();
 
            // Now we have the coeff size, we can determine the output size
            // Check it fits if necessary: symbol code,power, num,den of "coeff_size",
            // +1 for total coeff size
            uint64_t current_size = prev_size + ws + 1 + 2*ABS(coeff_size) + 1;
            if ( i > 0 ) {
                // and also symbols header, code,power of the symbol
                current_size += 4;
            }
            if ( write && must_fit_term && (sizeof(WORD)*current_size > (size_t)AM.MaxTer) ) {
                MLOCK(ErrorMessageLock);
                MesPrint("flint::to_argument_poly: output exceeds MaxTermSize");
                MUNLOCK(ErrorMessageLock);
                Terminate(-1);
            }
 
            WORD* term_size = 0;
            IFW(term_size = out++); ws++;
 
            if ( i > 0 ) {
                IFW(*out++ = SYMBOL); ws++;
                IFW(*out++ = 4); ws++; // The symbol array size, it is univariate
                IFW(*out++ = var_map_inv[0]); ws++;
                IFW(*out++ = i); ws++;
            }
 
            // Copy numerator
            for ( WORD j = 0; j < ABS(num_size); j++ ) {
                IFW(*out++ = tmp_coeff[j]); ws++;
            }
            for ( WORD j = ABS(num_size); j < ABS(coeff_size); j++ ) {
                IFW(*out++ = 0); ws++;
            }
            // Copy denominator
            for ( WORD j = 0; j < ABS(den_size); j++ ) {
                IFW(*out++ = tmp_den[j]); ws++;
            }
            for ( WORD j = ABS(den_size); j < ABS(coeff_size); j++ ) {
                IFW(*out++ = 0); ws++;
            }
 
            IFW(*out = 2*ABS(coeff_size) + 1); // the size of the coefficient
            IFW(if ( coeff_size < 0 ) { *out = -(*out); });
            IFW(out++); ws++;
 
            IFW(*term_size = out - term_size);
        }
 
    }
 
    if ( with_arghead ) {
        IFW(*arg_size = out - arg_size);
    }
    else {
        // with no arghead, we write a terminating zero
        IFW(*out++ = 0); ws++;
    }
 
    NumberFree(tmp_coeff, "flint::to_argument_poly");
    NumberFree(tmp_den, "flint::to_argument_poly");
 
    return ws;
}
 
// If no denscale argument is supplied, just set it to 1 and call the usual function
uint64_t flint::to_argument_poly(PHEAD WORD *out, const bool with_arghead,
    const bool must_fit_term, const bool write, const uint64_t prev_size, const fmpz_poly_t poly,
    const var_map_t &var_map) {
 
    flint::fmpz tmp;
    fmpz_set_ui(tmp.d, 1);
 
    uint64_t ret = flint::to_argument_poly(BHEAD out, with_arghead, must_fit_term, write, prev_size,
        poly, var_map, tmp.d);
 
    return ret;
}
/*
    #] flint::to_argument_poly :
*/
 
// Utility functions
/*
    #[ flint::util::simplify_fmpz :
*/
// Divide the GCD out of num and den
inline void flint::util::simplify_fmpz(fmpz_t num, fmpz_t den, fmpz_t gcd) {
    fmpz_gcd(gcd, num, den);
    if ( !fmpz_is_one(gcd) ) {
        fmpz_divexact(num, num, gcd);
        fmpz_divexact(den, den, gcd);
    }
}
/*
    #] flint::util::simplify_fmpz :
    #[ flint::util::simplify_fmpz_poly :
*/
// Divide the GCD out of num and den
inline void flint::util::simplify_fmpz_poly(fmpz_poly_t num, fmpz_poly_t den, fmpz_poly_t gcd) {
    fmpz_poly_gcd(gcd, num, den);
    if ( !fmpz_poly_is_one(gcd) ) {
#if __FLINT_RELEASE >= 30100
        // This should be faster than fmpz_poly_div, see https://github.com/flintlib/flint/pull/1766
        fmpz_poly_divexact(num, num, gcd);
        fmpz_poly_divexact(den, den, gcd);
#else
        fmpz_poly_div(num, num, gcd);
        fmpz_poly_div(den, den, gcd);
#endif
    }
}
/*
    #] flint::util::simplify_fmpz_poly :
    #[ flint::util::fix_sign_fmpz_mpoly_ratfun :
*/
inline void flint::util::fix_sign_fmpz_mpoly_ratfun(fmpz_mpoly_t num, fmpz_mpoly_t den,
    const fmpz_mpoly_ctx_t ctx) {
    // Fix sign to align with poly: the leading denominator term should have a positive coeff
    if ( fmpz_sgn(fmpz_mpoly_term_coeff_ref(den, 0, ctx)) == -1 ) {
        fmpz_mpoly_neg(num, num, ctx);
        fmpz_mpoly_neg(den, den, ctx);
    }
}
/*
    #] flint::util::fix_sign_fmpz_mpoly_ratfun :
    #[ flint::util::fix_sign_fmpz_poly_ratfun :
*/
inline void flint::util::fix_sign_fmpz_poly_ratfun(fmpz_poly_t num, fmpz_poly_t den) {
    // Fix sign to align with poly: the leading denominator term should have a positive coeff
    if ( fmpz_sgn(fmpz_poly_get_coeff_ptr(den, fmpz_poly_degree(den))) == -1 ) {
        fmpz_poly_neg(num, num);
        fmpz_poly_neg(den, den);
    }
}
/*
    #] flint::util::fix_sign_fmpz_poly_ratfun :
*/