poly.h

#pragma once
/* #[ 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 : */ 
 
extern "C" {
#include "form3.h"
}
 
#include <iostream>
#include <string>
#include <vector>
 
// for debugging
// #define POLY_MOVE_DEBUG
 
// macros for tform
#ifndef WITHPTHREADS
#define POLY_GETIDENTITY(X)
#define POLY_STOREIDENTITY
#else
#define POLY_GETIDENTITY(X) ALLPRIVATES *B = (X).Bpointer
#define POLY_STOREIDENTITY Bpointer = B
#endif
 
// maximum size of the hash table used for multiplication and division
 
const int POLY_MAX_HASH_SIZE = MiN(1<<20, MAXPOSITIVE);
 
class poly {
 
public:
 
    // variables
#ifdef WITHPTHREADS
    ALLPRIVATES *Bpointer;
#endif
 
    WORD *terms;
    LONG size_of_terms;
    WORD modp,modn;
 
    // constructors/destructor
    poly (PHEAD int, WORD=-1, WORD=1);
    poly (PHEAD const UWORD *, WORD, WORD=-1, WORD=1);
    poly (const poly &, WORD=-1, WORD=1);
    ~poly ();
 
    // operators
    poly& operator+= (const poly&);
    poly& operator-= (const poly&);
    poly& operator*= (const poly&);
    poly& operator/= (const poly&);
    poly& operator%= (const poly&);
 
    const poly operator+ (const poly&) const;
    const poly operator- (const poly&) const;
    const poly operator* (const poly&) const;
    const poly operator/ (const poly&) const;
    const poly operator% (const poly&) const;
 
    bool operator== (const poly&) const;
    bool operator!= (const poly&) const;
    poly& operator= (const poly &);
    WORD& operator[] (int);
    const WORD& operator[] (int) const;
 
#if __cplusplus >= 201103L
    // support move semantics
    poly (poly &&, WORD=-1, WORD=1) noexcept;
    poly& operator= (poly &&) noexcept;
#endif
 
    // memory management
    void termscopy (const WORD *, int, int);
    void check_memory(int);
    void expand_memory(int);
 
    // check type of polynomial
    bool is_zero () const;
    bool is_one () const;
    bool is_integer () const;
    bool is_monomial () const;
    int is_dense_univariate () const;
 
    // properties
    int sign () const;
    int degree (int) const;
    int total_degree () const;
    int first_variable () const;
    int number_of_terms () const;
    const std::vector<int> all_variables () const;  
    const poly integer_lcoeff () const;
    const poly lcoeff_univar (int) const;
    const poly lcoeff_multivar (int) const;
    const poly coefficient (int, int) const;
    const poly derivative (int) const;
    
    // modulo calculus
    void setmod(WORD, WORD=1);
    void coefficients_modulo (UWORD *, WORD, bool);
 
    // simple polynomials
    static const poly simple_poly (PHEAD int, int=0, int=1, int=0, int=1);
    static const poly simple_poly (PHEAD int, const poly&, int=1, int=0, int=1);
 
    // conversion from/to form notation
    static void get_variables (PHEAD const std::vector<WORD *> &, bool, bool);
    static const poly argument_to_poly (PHEAD WORD *, bool, bool, poly *den=NULL);
    static void poly_to_argument (const poly &, WORD *, bool);
    static void poly_to_argument_with_den (const poly &, WORD, const UWORD *, WORD *, bool);
    int size_of_form_notation () const;
    int size_of_form_notation_with_den (WORD) const;
    const poly & normalize ();
 
    // operations for coefficient lists
    static const poly from_coefficient_list (PHEAD const std::vector<WORD> &, int, WORD);
    static const std::vector<WORD> to_coefficient_list (const poly &);
    static const std::vector<WORD> coefficient_list_divmod (const std::vector<WORD> &, const std::vector<WORD> &, WORD, int);
 
    // string output for debugging
    const std::string to_string() const;
 
    // monomials
    static int monomial_compare (PHEAD const WORD *, const WORD *);
    WORD last_monomial_index () const;
    WORD* last_monomial () const;
    int compare_degree_vector(const poly &) const;
    std::vector<int> degree_vector() const;
    int compare_degree_vector(const std::vector<int> &) const;
 
    // (internal) mathematical operations
    static void add (const poly &, const poly &, poly &);
    static void sub (const poly &, const poly &, poly &);
    static void mul (const poly &, const poly &, poly &);
    static void div (const poly &, const poly &, poly &);
    static void mod (const poly &, const poly &, poly &);
    static void divmod (const poly &, const poly &, poly &, poly &, bool only_divides);
    static bool divides (const poly &, const poly &);
    
    static void mul_one_term (const poly &, const poly &, poly &);
    static void mul_univar (const poly &, const poly &, poly &, int);
    static void mul_heap (const poly &, const poly &, poly &);
 
    static void divmod_one_term (const poly &, const poly &, poly &, poly &, bool);
    static void divmod_univar (const poly &, const poly &, poly &, poly &, int, bool);
    static void divmod_heap (const poly &, const poly &, poly &, poly &, bool, bool, bool&);
    
    static void push_heap (PHEAD WORD **, int);
    static void pop_heap (PHEAD WORD **, int);
};
 
// comparison class for monomials (for std::sort)
class monomial_larger {
 
public:
    
#ifndef WITHPTHREADS
    monomial_larger() {}
#else
    ALLPRIVATES *B;
    monomial_larger(ALLPRIVATES *b): B(b) {}
#endif
    
    bool operator()(const WORD *a, const WORD *b) {
        return poly::monomial_compare(BHEAD a, b) > 0;
    }
};
 
// stream operator
std::ostream& operator<< (std::ostream &, const poly &);
 
// inline function definitions
 
#if __cplusplus >= 201103L
 
inline poly::poly (poly &&a, WORD modp, WORD modn) noexcept {
#ifdef POLY_MOVE_DEBUG
    std::cout << "poly move ctor" << std::endl;
#endif
    POLY_GETIDENTITY(a);
    POLY_STOREIDENTITY;
    terms = a.terms;
    size_of_terms = a.size_of_terms;
    this->modp = a.modp;
    this->modn = a.modn;
    if (modp != -1) {
        setmod(modp,modn);
    }
    a.terms = nullptr;
    a.size_of_terms = -1;
}
 
inline poly& poly::operator= (poly &&a) noexcept {
#ifdef POLY_MOVE_DEBUG
    std::cout << "poly move assign" << std::endl;
#endif
    if (this != &a) {
        WORD *old_terms = terms;
        LONG old_size_of_terms = size_of_terms;
        terms = a.terms;
        size_of_terms = a.size_of_terms;
        modp = a.modp;
        modn = a.modn;
        a.terms = old_terms;
        a.size_of_terms = old_size_of_terms;
    }
    return *this;
}
 
#endif
 
/*   Checks whether the terms array is large enough to add another
 *   term (of size AM.MaxTal) to the polynomials. In case not, it is
 *   expanded.
 */
inline void poly::check_memory (int i) {
    POLY_GETIDENTITY(*this);
    if (i + 3 + AN.poly_num_vars + AM.MaxTal >= size_of_terms) expand_memory(i + AM.MaxTal);
//  Used to be i+2 but there should also be space for a trailing zero
}
 
// indexing operators
inline WORD& poly::operator[] (int i) {
    return terms[i];
}
 
inline const WORD& poly::operator[] (int i) const {
    return terms[i];
}
 
/*   Copies "num" WORD-sized terms from the pointer "source" to the
 *   current polynomial at index "dest"
 */
inline void poly::termscopy (const WORD *source, int dest, int num) {
#ifdef POLY_MOVE_DEBUG
    std::cout << "poly termscopy: num=" << num << std::endl;
#endif
    memcpy (terms+dest, source, num*sizeof(WORD));
}