src/AlgebraicCore/PPMonoidEv.C

//   Copyright (c)  2004-2011  John Abbott

//   This file is part of the source of CoCoALib, the CoCoA Library.

//   CoCoALib 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.

//   CoCoALib 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 CoCoALib.  If not, see <http://www.gnu.org/licenses/>.


// Implementation of class PPMonoidEvImpl

#include "CoCoA/PPMonoidEv.H"

#include "CoCoA/BigIntOps.H"
#include "CoCoA/DivMask.H"
#include "CoCoA/MemPool.H"
#include "CoCoA/PPMonoid.H"
#include "CoCoA/PPOrdering.H"
#include "CoCoA/assert.H"
#include "CoCoA/convert.H"
#include "CoCoA/degree.H"
#include "CoCoA/error.H"
#include "CoCoA/matrix.H"
#include "CoCoA/symbol.H"

#include <algorithm>
using std::min;
using std::max;
//using std::swap;
#include <cstring>
using std::memcpy;
#include <iostream>
using std::ostream;
#include<limits>
using std::numeric_limits;
#include <memory>
using std::unique_ptr;
#include <vector>
using std::vector;


namespace CoCoA
{

  /*-- class PPMonoidEvSmallExpImpl ---------------------------------------*/
  /*-- class PPMonoidEvImpl ---------------------------------------*/
  /**

  \brief Implementation power product monoid with exponent vector

  PPMonoidEvSmallExpImpl implements a power product monoid for safe testing.
  It stores the exponents as a C vector of SmallExponent_t.
  It is not designed to be efficient (but it might be ;-)

  */
  /*-----------------------------------------------------------------*/

  class PPMonoidEvSmallExpImpl: public PPMonoidBase
  {
    typedef PPMonoidElemRawPtr RawPtr;           // just to save typing
    typedef PPMonoidElemConstRawPtr ConstRawPtr; // just to save typing

    static const unsigned long ourMaxExp;        // defined below; value is just numeric_limits<SmallExponent_t>::max()

    class CmpBase
    {
    public:
      CmpBase(long NumIndets, long GradingDim);
      virtual ~CmpBase();
    public:
      virtual int myCmpExpvs(const SmallExponent_t* v1, const SmallExponent_t* v2) const =0;
      virtual void myWDeg(degree& d, const SmallExponent_t* v) const =0;
      //      virtual int myCmpWDeg(const SmallExponent_t* v1, const SmallExponent_t* v2) const =0;
      virtual int myCmpWDegPartial(const SmallExponent_t* v1, const SmallExponent_t* v2, long) const =0;
      int myCmpWDeg(const SmallExponent_t* v1, const SmallExponent_t* v2) const { return myCmpWDegPartial(v1, v2, myGradingDim); }

    protected:
      long myNumIndets;
      long myGradingDim;
    };

    class LexImpl: public CmpBase
    {
    public:
      LexImpl(long NumIndets);
      virtual int myCmpExpvs(const SmallExponent_t* v1, const SmallExponent_t* v2) const;
      virtual void myWDeg(degree& d, const SmallExponent_t* v) const;
      virtual int myCmpWDegPartial(const SmallExponent_t* v1, const SmallExponent_t* v2, long) const;
    };


    class StdDegLexImpl: public CmpBase
    {
    public:
      StdDegLexImpl(long NumIndets);
      virtual int myCmpExpvs(const SmallExponent_t* v1, const SmallExponent_t* v2) const;
      virtual void myWDeg(degree& d, const SmallExponent_t* v) const;
      virtual int myCmpWDegPartial(const SmallExponent_t* v1, const SmallExponent_t* v2, long) const;
    };


    class StdDegRevLexImpl: public CmpBase
    {
    public:
      StdDegRevLexImpl(long NumIndets);
      virtual int myCmpExpvs(const SmallExponent_t* v1, const SmallExponent_t* v2) const;
      virtual void myWDeg(degree& d, const SmallExponent_t* v) const;
      virtual int myCmpWDegPartial(const SmallExponent_t* v1, const SmallExponent_t* v2, long) const;
    };


    class MatrixOrderingImpl: public CmpBase
    {
    public:
      MatrixOrderingImpl(long NumIndets, long GradingDim, const ConstMatrixView& OrderMatrix);
      virtual int myCmpExpvs(const SmallExponent_t* v1, const SmallExponent_t* v2) const;
      virtual void myWDeg(degree& d, const SmallExponent_t* v) const;
      virtual int myCmpWDegPartial(const SmallExponent_t* v1, const SmallExponent_t* v2, long) const;
    private:
      std::vector< std::vector<BigInt> > myOrderMatrix;
    };


  public:
    PPMonoidEvSmallExpImpl(const std::vector<symbol>& IndetNames, const PPOrdering& ord);
    ~PPMonoidEvSmallExpImpl();
  private: // disable copy ctor and assignment
    explicit PPMonoidEvSmallExpImpl(const PPMonoidEvSmallExpImpl& copy);  // NEVER DEFINED -- copy ctor disabled
    PPMonoidEvSmallExpImpl& operator=(const PPMonoidEvSmallExpImpl& rhs); // NEVER DEFINED -- assignment disabled

  public:
    void contents() const; // FOR DEBUGGING ONLY

    const std::vector<PPMonoidElem>& myIndets() const;               ///< std::vector whose n-th entry is n-th indet as PPMonoidElem

    // The functions below are operations on power products owned by PPMonoidEvSmallExpImpl
    const PPMonoidElem& myOne() const;
    using PPMonoidBase::myNew;    // disable warnings of overloading
    PPMonoidElemRawPtr myNew() const;                                ///< ctor from nothing
    PPMonoidElemRawPtr myNew(PPMonoidElemConstRawPtr rawpp) const;   ///< ctor by assuming ownership
    PPMonoidElemRawPtr myNew(const std::vector<long>& expv) const;   ///< ctor from exp vector
    void myDelete(RawPtr rawpp) const;                               ///< dtor, frees pp
    void mySwap(RawPtr rawpp1, RawPtr rawpp2) const;                 ///< swap(pp1, pp2);
    void myAssignOne(RawPtr rawpp) const;                            ///< pp = 1
    void myAssign(RawPtr rawpp, ConstRawPtr rawpp1) const;           ///< p = pp1
    void myAssign(RawPtr rawpp, const std::vector<long>& expv) const;///< pp = expv (assign from exp vector, all exps >= 0)

    void myMul(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const;   ///< pp = pp1*pp2
    using PPMonoidBase::myMulIndetPower;    // disable warnings of overloading
    void myMulIndetPower(RawPtr rawpp, long indet, long exp) const;           ///< pp *= indet^exp, assumes exp >= 0
    void myDiv(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const;   ///< pp = pp1/pp2
    void myColon(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const; ///< pp = pp1/gcd(pp1,pp2)
    void myGcd(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const;   ///< pp = gcd(pp1,pp2)
    void myLcm(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const;   ///< pp = lcm(pp1,pp2)
    void myRadical(RawPtr rawpp, ConstRawPtr rawpp1) const;                   ///< pp = radical(pp1)
    void myPowerSmallExp(RawPtr rawpp, ConstRawPtr rawpp1, long exp) const;   ///< pp = pp1^exp, assumes exp >= 0
    void myPowerOverflowCheck(ConstRawPtr rawpp1, long exp) const;            ///< throw if pp1^exp would overflow, assumes exp >= 0


    bool myIsOne(ConstRawPtr rawpp) const;                                   ///< is pp = 1?
    bool myIsIndet(long& index, ConstRawPtr rawpp) const;                    ///< true iff pp is an indet
    bool myIsCoprime(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const;          ///< are pp1 & pp2 coprime?
    bool myIsEqual(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const;            ///< is pp1 equal to pp2?
    bool myIsDivisible(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const;        ///< does pp2 divide pp1?
    bool myIsSqFree(ConstRawPtr rawpp) const;                                ///< is pp equal to its radical?

    int myCmp(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const;                 ///< -1,0,1 as pp1 < = > pp2
    long myStdDeg(ConstRawPtr rawpp) const;                                  ///< standard degree of pp
    void myWDeg(degree& d, ConstRawPtr rawpp) const;                         ///< d = grading(pp)
    int myCmpWDeg(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const;             ///< <0, =0, >0 as wdeg(pp1) < = > wdeg(pp2)
    int myCmpWDegPartial(ConstRawPtr rawpp1, ConstRawPtr rawpp2, long) const; ///< as myCmpWDeg wrt the first weights
    long myExponent(ConstRawPtr rawpp, long indet) const;                    ///< exponent of indet in pp
    void myBigExponent(BigInt& EXP, ConstRawPtr rawpp, long indet) const;    ///< EXP = exponent of indet in pp
    void myExponents(std::vector<long>& expv, ConstRawPtr rawpp) const;      ///< expv[i] = exponent(pp,i)
    void myBigExponents(std::vector<BigInt>& v, ConstRawPtr rawpp) const;    ///< get exponents, SHOULD BE myExponents ???
    void myIndetsIn(std::vector<bool>& v, ConstRawPtr rawpp) const;          ///< v[i] = true if i-th indet has exponent != 0
    void myOutputSelf(std::ostream& out) const;                              ///< out << PPM
    // INHERITED DEFINITION of virtual  void myOutput(std::ostream& out, ConstRawPtr rawpp) const;
    void myDebugPrint(std::ostream& out, ConstRawPtr rawpp) const;           ///< print pp in debugging format???


  private: // auxiliary functions
    SmallExponent_t* myExpv(RawPtr) const;
    const SmallExponent_t* myExpv(ConstRawPtr) const;

    void myComputeDivMask(DivMask& dm, const DivMaskRule& DivMaskImpl, ConstRawPtr rawpp) const; ///< used by PPWithMask
    bool myCheckExponents(const std::vector<long>& expv) const;

  private: // data members
    ///@name Class members
    //@{
    const long myEntrySize;       ///< size in ?bytes???
    mutable MemPool myMemMgr;     // IMPORTANT: this must come *before* myIndetVector and myOnePtr.
    vector<PPMonoidElem> myIndetVector; ///< the indets as PPMonoidElems
    unique_ptr<CmpBase> myOrdPtr;         ///< actual implementation of the ordering [should be const???]
    unique_ptr<PPMonoidElem> myOnePtr;
    //@}
  };

  // static variable
  const unsigned long PPMonoidEvSmallExpImpl::ourMaxExp = numeric_limits<SmallExponent_t>::max();


  // dtor not inline for some compiler (?) on PPC
  PPMonoidEvSmallExpImpl::CmpBase::~CmpBase() {}


  // File local inline functions

  inline SmallExponent_t* PPMonoidEvSmallExpImpl::myExpv(RawPtr rawpp) const
  {
    return static_cast<SmallExponent_t*>(rawpp.myRawPtr());
  }


  inline const SmallExponent_t* PPMonoidEvSmallExpImpl::myExpv(ConstRawPtr rawpp) const
  {
    return static_cast<const SmallExponent_t*>(rawpp.myRawPtr());
  }


  bool PPMonoidEvSmallExpImpl::myCheckExponents(const std::vector<long>& expv) const
  {
    // Check expv.size == myNumIndets.
    // Check exps are non-neg and not too big.
    if (len(expv) != myNumIndets) return false;
    for (long i=0; i < myNumIndets; ++i)
      if (expv[i] < 0 || static_cast<unsigned long>(expv[i]) > numeric_limits<SmallExponent_t>::max()) return false;
    return true;
  }


  //----   Constructors & destructor   ----//

  PPMonoidEvSmallExpImpl::PPMonoidEvSmallExpImpl(const std::vector<symbol>& IndetNames, const PPOrdering& ord):
      PPMonoidBase(ord, IndetNames),
      myEntrySize(sizeof(SmallExponent_t)*myNumIndets),
      myMemMgr(myEntrySize, "PPMonoidEvSmallExpImpl.myMemMgr"),
      myIndetVector()
  {
    // Put the correct implementation in myOrdPtr...  [VERY UGLY CODE!!!]
    if (IsLex(ord))
      myOrdPtr.reset(new LexImpl(myNumIndets));
    else if (IsStdDegLex(ord))
      myOrdPtr.reset(new StdDegLexImpl(myNumIndets));
    else if (IsStdDegRevLex(ord))
      myOrdPtr.reset(new StdDegRevLexImpl(myNumIndets));
    else
      myOrdPtr.reset(new MatrixOrderingImpl(myNumIndets, GradingDim(ord), OrdMat(ord)));

    myRefCountInc();  // this is needed for exception cleanliness, in case one of the lines below throws
    myOnePtr.reset(new PPMonoidElem(PPMonoid(this)));
    {
      // IMPORTANT: this block destroys pp *before* the call to myRefCountZero.
      PPMonoidElem pp(PPMonoid(this));
      vector<long> expv(myNumIndets);
      for (long i=0; i < myNumIndets; ++i)
      {
        expv[i] = 1;
        myAssign(raw(pp), expv);
        myIndetVector.push_back(pp);
        expv[i] = 0;
      }
    }
    myRefCountZero();
  }


  PPMonoidEvSmallExpImpl::~PPMonoidEvSmallExpImpl()
  {}

/////////////////////////////////////////////////////////////////////////////


  const std::vector<PPMonoidElem>& PPMonoidEvSmallExpImpl::myIndets() const
  {
    return myIndetVector;
  }


  const PPMonoidElem& PPMonoidEvSmallExpImpl::myOne() const
  {
    return *myOnePtr;
  }


  PPMonoidElemRawPtr PPMonoidEvSmallExpImpl::myNew() const
  {
    PPMonoidElemRawPtr rawpp(myMemMgr.alloc());
    myAssignOne(rawpp); // cannot throw
    return rawpp;
  }

  PPMonoidElemRawPtr PPMonoidEvSmallExpImpl::myNew(PPMonoidElemConstRawPtr rawcopypp) const
  {
    PPMonoidElemRawPtr rawpp(myMemMgr.alloc());
    myAssign(rawpp, rawcopypp); // cannot throw
    return rawpp;
  }


  PPMonoidElemRawPtr PPMonoidEvSmallExpImpl::myNew(const std::vector<long>& expv) const
  {
    CoCoA_ASSERT(myCheckExponents(expv));
    PPMonoidElemRawPtr rawpp(myMemMgr.alloc());
    myAssign(rawpp, expv); // cannot throw
    return rawpp;
  }


  void PPMonoidEvSmallExpImpl::myAssignOne(RawPtr rawpp) const
  {
    SmallExponent_t* const expv = myExpv(rawpp);
    for (long i = 0; i < myNumIndets; ++i)
      expv[i] = 0;
  }


  void PPMonoidEvSmallExpImpl::myAssign(RawPtr rawpp, ConstRawPtr rawpp1) const
  {
    if (rawpp == rawpp1) return;

//     SmallExponent_t* const expv = myExpv(rawpp);
//     const SmallExponent_t* const expv1 = myExpv(rawpp1);
//     std::copy(expv1, expv1+myNumIndets, expv);
    memcpy(myExpv(rawpp), myExpv(rawpp1), myNumIndets*sizeof(SmallExponent_t));
  }

  void PPMonoidEvSmallExpImpl::myAssign(RawPtr rawpp, const vector<long>& expv1) const
  {
    CoCoA_ASSERT(myCheckExponents(expv1));

    SmallExponent_t* const expv = myExpv(rawpp);
    for (long i = 0; i < myNumIndets; ++i)
      expv[i] = NumericCast<SmallExponent_t>(expv1[i]);
  }


  void PPMonoidEvSmallExpImpl::myDelete(RawPtr rawpp) const
  {
    myMemMgr.free(rawpp.myRawPtr());
  }


  void PPMonoidEvSmallExpImpl::mySwap(RawPtr rawpp1, RawPtr rawpp2) const
  {
    if (rawpp1 == rawpp2) return;
    SmallExponent_t* const expv1 = myExpv(rawpp1);
    SmallExponent_t* const expv2 = myExpv(rawpp2);
    for (long i = 0; i < myNumIndets; ++i)
      std::swap(expv1[i], expv2[i]);
  }


  void PPMonoidEvSmallExpImpl::myMul(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    // No worries about aliasing.
    SmallExponent_t* const expv = myExpv(rawpp);
    const SmallExponent_t* const expv1 = myExpv(rawpp1);
    const SmallExponent_t* const expv2 = myExpv(rawpp2);
    for (long i=0; i < myNumIndets; ++i)
    {
      CoCoA_ASSERT("Exponent Overflow" && expv1[i] <= std::numeric_limits<SmallExponent_t>::max()-expv2[i]);
      expv[i] = expv1[i] + expv2[i];
    }
  }


  void PPMonoidEvSmallExpImpl::myMulIndetPower(RawPtr rawpp, long indet, long exp) const  // assumes exp >= 0
  {
    CoCoA_ASSERT(exp >= 0);
    CoCoA_ASSERT(0 <= indet && indet < myNumIndets);
    SmallExponent_t* const expv = myExpv(rawpp);
    // If CoCoA_DEBUG active, check for exponent overflow...
    CoCoA_ASSERT("Exponent Overflow" && ourMaxExp - expv[indet] >= static_cast<unsigned long>(exp));
    expv[indet] += static_cast<SmallExponent_t>(exp);  // cast to keep M$ compiler quiet
  }


  void PPMonoidEvSmallExpImpl::myDiv(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    // No worries about aliasing.
    SmallExponent_t* const expv = myExpv(rawpp);
    const SmallExponent_t* const expv1 = myExpv(rawpp1);
    const SmallExponent_t* const expv2 = myExpv(rawpp2);
    for (long i=0; i < myNumIndets; ++i)
    {
      CoCoA_ASSERT("Exponent Underflow" && expv1[i] >= expv2[i]);
      expv[i] = expv1[i] - expv2[i];
    }
  }


  void PPMonoidEvSmallExpImpl::myColon(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    // No worries about aliasing.
    SmallExponent_t* const expv = myExpv(rawpp);
    const SmallExponent_t* const expv1 = myExpv(rawpp1);
    const SmallExponent_t* const expv2 = myExpv(rawpp2);

    for (long i = 0; i < myNumIndets; ++i)
      if (expv1[i] > expv2[i])
        expv[i] = expv1[i] - expv2[i];
      else
        expv[i] = 0;
  }


  void PPMonoidEvSmallExpImpl::myGcd(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    // No worries about aliasing.
    SmallExponent_t* const expv = myExpv(rawpp);
    const SmallExponent_t* const expv1 = myExpv(rawpp1);
    const SmallExponent_t* const expv2 = myExpv(rawpp2);

    for (long i = 0; i < myNumIndets; ++i)
      expv[i] = min(expv1[i], expv2[i]);
  }


  void PPMonoidEvSmallExpImpl::myLcm(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    // No worries about aliasing.
    SmallExponent_t* const expv = myExpv(rawpp);
    const SmallExponent_t* const expv1 = myExpv(rawpp1);
    const SmallExponent_t* const expv2 = myExpv(rawpp2);

    for (long i = 0; i < myNumIndets; ++i)
      expv[i] = max(expv1[i], expv2[i]);
  }


  void PPMonoidEvSmallExpImpl::myRadical(RawPtr rawpp, ConstRawPtr rawpp1) const
  {
    SmallExponent_t* const expv = myExpv(rawpp);
    const SmallExponent_t* const expv1 = myExpv(rawpp1);

    for (long i = 0; i < myNumIndets; ++i)
      expv[i] = (expv1[i] > 0);
  }


  void PPMonoidEvSmallExpImpl::myPowerSmallExp(RawPtr rawpp, ConstRawPtr rawpp1, long LongExp) const  // assumes exp >= 0
  {
    CoCoA_ASSERT(LongExp >= 0);
#ifdef CoCoA_DEBUG
    myPowerOverflowCheck(rawpp1, LongExp);
#endif
    if (static_cast<unsigned long>(LongExp) > ourMaxExp)
      CoCoA_ERROR(ERR::ExpTooBig, "PPMonoidEvSmallExpImpl::myPowerSmallExp");
    const SmallExponent_t exp = static_cast<SmallExponent_t>(LongExp);

    SmallExponent_t* const expv = myExpv(rawpp);
    const SmallExponent_t* const expv1 = myExpv(rawpp1);
    for (long i = 0; i < myNumIndets; ++i)
      expv[i] = exp * expv1[i];
  }


  void PPMonoidEvSmallExpImpl::myPowerOverflowCheck(ConstRawPtr rawpp, long LongExp) const
  {
    if (LongExp == 0 || LongExp == 1) return;
    CoCoA_ASSERT(LongExp >= 0);
    const char* const FnName = "PPMonoidEvSmallExpImpl::myPowerOverflowCheck";
    if (static_cast<unsigned long>(LongExp) > ourMaxExp)
      CoCoA_ERROR(ERR::ExpTooBig, FnName);
    const SmallExponent_t exp = static_cast<SmallExponent_t>(LongExp);
    const SmallExponent_t limit = ourMaxExp/exp;

    const SmallExponent_t* const expv = myExpv(rawpp);
    for (long i = 0; i < myNumIndets; ++i)
    {
      if (expv[i] > limit)
        CoCoA_ERROR(ERR::ExpTooBig, FnName);
    }
  }


  bool PPMonoidEvSmallExpImpl::myIsOne(ConstRawPtr rawpp) const
  {
    const SmallExponent_t* const expv = myExpv(rawpp);

    for (long i = 0; i < myNumIndets; ++i)
      if (expv[i] != 0) return false;

    return true;
  }


  bool PPMonoidEvSmallExpImpl::myIsIndet(long& index, ConstRawPtr rawpp) const
  {
    const SmallExponent_t* const expv = myExpv(rawpp);
    long j = myNumIndets;
    for (long i = 0; i < myNumIndets; ++i)
    {
      if (expv[i] == 0) continue;
      if (j != myNumIndets || expv[i] != 1) return false;
      j = i;
    }
    if (j == myNumIndets) return false;
    index = j;
    return true;
  }


  bool PPMonoidEvSmallExpImpl::myIsCoprime(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    const SmallExponent_t* const expv1 = myExpv(rawpp1);
    const SmallExponent_t* const expv2 = myExpv(rawpp2);

    for (long i = 0; i < myNumIndets; ++i)
      if (expv1[i] != 0 && expv2[i] != 0) return false;

    return true;
  }


  bool PPMonoidEvSmallExpImpl::myIsEqual(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    const SmallExponent_t* const expv1 = myExpv(rawpp1);
    const SmallExponent_t* const expv2 = myExpv(rawpp2);

    for (long i = 0; i < myNumIndets; ++i)
      if (expv1[i] != expv2[i]) return false;
    return true;
  }


  bool PPMonoidEvSmallExpImpl::myIsDivisible(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    const SmallExponent_t* const expv1 = myExpv(rawpp1);
    const SmallExponent_t* const expv2 = myExpv(rawpp2);

    for (long i = 0; i < myNumIndets; ++i)
      if (expv1[i] < expv2[i]) return false;

    return true;
  }


  bool PPMonoidEvSmallExpImpl::myIsSqFree(ConstRawPtr rawpp) const
  {
    const SmallExponent_t* const expv = myExpv(rawpp);

    for (long i = 0; i < myNumIndets; ++i)
      if (expv[i] > 1) return false;

    return true;
  }


  int PPMonoidEvSmallExpImpl::myCmp(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    return myOrdPtr->myCmpExpvs(myExpv(rawpp1), myExpv(rawpp2));
  }


  long PPMonoidEvSmallExpImpl::myStdDeg(ConstRawPtr rawpp) const
  {
    const SmallExponent_t* const expv = myExpv(rawpp);
    long d=0;
    for (long i=0; i < myNumIndets; ++i)
      d += expv[i];
    return d;
  }


  void PPMonoidEvSmallExpImpl::myWDeg(degree& d, ConstRawPtr rawpp) const
  {
    myOrdPtr->myWDeg(d, myExpv(rawpp));
  }


  int PPMonoidEvSmallExpImpl::myCmpWDeg(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    return myOrdPtr->myCmpWDeg(myExpv(rawpp1), myExpv(rawpp2));
  }


  int PPMonoidEvSmallExpImpl::myCmpWDegPartial(ConstRawPtr rawpp1, ConstRawPtr rawpp2, long i) const
  {
    return myOrdPtr->myCmpWDegPartial(myExpv(rawpp1), myExpv(rawpp2), i);
  }


  long PPMonoidEvSmallExpImpl::myExponent(ConstRawPtr rawpp, long indet) const
  {
    CoCoA_ASSERT(indet < myNumIndets);
    return NumericCast<long>(myExpv(rawpp)[indet]);
  }

  void PPMonoidEvSmallExpImpl::myBigExponent(BigInt& EXP, ConstRawPtr rawpp, long indet) const
  {
    CoCoA_ASSERT(indet < myNumIndets);
    EXP = myExpv(rawpp)[indet];
  }


  void PPMonoidEvSmallExpImpl::myExponents(std::vector<long>& v, ConstRawPtr rawpp) const
  {
    CoCoA_ASSERT(len(v) == myNumIndets);
    const SmallExponent_t* const expv = myExpv(rawpp);
    for (long i=0; i < myNumIndets; ++i)
      v[i] = NumericCast<long>(expv[i]);
  }


  void PPMonoidEvSmallExpImpl::myBigExponents(std::vector<BigInt>& expv, ConstRawPtr rawpp) const
  {
    CoCoA_ASSERT(len(expv) == myNumIndets);
    const SmallExponent_t* const v = myExpv(rawpp);
    for (long i=0; i < myNumIndets; ++i)  expv[i] = v[i];
  }


  void PPMonoidEvSmallExpImpl::myIndetsIn(std::vector<bool>& v, ConstRawPtr rawpp) const
  {
    CoCoA_ASSERT(len(v) == myNumIndets);
    const SmallExponent_t* const expv = myExpv(rawpp);
    for (long i=0; i < myNumIndets; ++i)
      if (expv[i] != 0)
        v[i] = true;
  }

  void PPMonoidEvSmallExpImpl::myComputeDivMask(DivMask& dm, const DivMaskRule& DivMaskImpl, ConstRawPtr rawpp) const
  {
    DivMaskImpl->myAssignFromExpv(dm, myExpv(rawpp), myNumIndets);
  }


  void PPMonoidEvSmallExpImpl::myOutputSelf(std::ostream& out) const
  {
    if (!out) return;  // short-cut for bad ostreams
    out << "PPMonoidEv(" << myNumIndets << ", " << myOrd <<")";
  }


  void PPMonoidEvSmallExpImpl::myDebugPrint(std::ostream& out, ConstRawPtr rawpp) const
  {
    if (!out) return;  // short-cut for bad ostreams
    out << "DEBUG PP: myNumIndets=" << myNumIndets << ", exps=[";
    for (long i=0; i < myNumIndets; ++i)
      out << myExponent(rawpp, i) << " ";
    out << "]" << std::endl;
  }


//--- orderings ----------------------------------------

  PPMonoidEvSmallExpImpl::CmpBase::CmpBase(long NumIndets, long GradingDim):
    myNumIndets(NumIndets),
    myGradingDim(GradingDim)
  {
    CoCoA_ASSERT(NumIndets > 0);
    CoCoA_ASSERT(NumIndets < 1000000); // complain about ridiculously large number of indets
    CoCoA_ASSERT(GradingDim >= 0);
    CoCoA_ASSERT(GradingDim <= NumIndets);
  }


//--- LexImpl

  PPMonoidEvSmallExpImpl::LexImpl::LexImpl(long NumIndets):
    CmpBase(NumIndets, 0)
  {
  }


  int PPMonoidEvSmallExpImpl::LexImpl::myCmpExpvs(const SmallExponent_t* v1, const SmallExponent_t* v2) const
  {
    for (long i=0; i<myNumIndets; ++i)
      if (v1[i] != v2[i])
      {
        if (v1[i] > v2[i]) return 1; else return -1;
      }
    return 0;
  }


  void PPMonoidEvSmallExpImpl::LexImpl::myWDeg(degree& /*d*/, const SmallExponent_t* /*v*/) const
  {
    // deliberately does nothing because GradingDim=0
    //??? should it assign:  d = []
  }


  int PPMonoidEvSmallExpImpl::LexImpl::myCmpWDegPartial(const SmallExponent_t* /*v1*/, const SmallExponent_t* /*v2*/, long /*i*/) const
  {
    return 0; // GradingDim=0, and all degrees in Z^0 are equal
  }

//--- StdDegLexImpl

  PPMonoidEvSmallExpImpl::StdDegLexImpl::StdDegLexImpl(long NumIndets):
    CmpBase(NumIndets, 1)
  {
  }


  void PPMonoidEvSmallExpImpl::StdDegLexImpl::myWDeg(degree& d, const SmallExponent_t* v) const
  {
    long deg=0;
    for (long i=0; i < myNumIndets; ++i)
      deg += v[i];
    SetComponent(d, 0, deg);
  }


  int PPMonoidEvSmallExpImpl::StdDegLexImpl::myCmpExpvs(const SmallExponent_t* v1, const SmallExponent_t* v2) const
  {
    long deg1=0, deg2=0;
    for (long i=0; i < myNumIndets; ++i)
    {
      deg1 += v1[i];
      deg2 += v2[i];
    }
    if (deg1 != deg2)
    {
      if (deg1 > deg2) return 1; else return -1;
    }

    for (long i=0; i < myNumIndets; ++i)
      if (v1[i] != v2[i])
      {
        if (v1[i] > v2[i]) return 1; else return -1;
      }
    return 0;
  }


  int PPMonoidEvSmallExpImpl::StdDegLexImpl::myCmpWDegPartial(const SmallExponent_t* v1, const SmallExponent_t* v2, long i) const
  {
    if (i==0) return 0; // ie: wrt to 0-grading
    long deg1=0, deg2=0;
    for (long i=0; i < myNumIndets; ++i)
    {
      deg1 += v1[i];
      deg2 += v2[i];
    }
    if (deg1 != deg2)
    {
      if (deg1 > deg2) return 1; else return -1;
    }
    return 0;
  }


//--- StdDegRevLexImpl

  PPMonoidEvSmallExpImpl::StdDegRevLexImpl::StdDegRevLexImpl(long NumIndets):
    CmpBase(NumIndets, 1)
  {
  }


  void PPMonoidEvSmallExpImpl::StdDegRevLexImpl::myWDeg(degree& d, const SmallExponent_t* v) const
  {
    long deg=0;
    for (long i=0; i < myNumIndets; ++i)
      deg += v[i];
    SetComponent(d, 0, deg);
  }


  int PPMonoidEvSmallExpImpl::StdDegRevLexImpl::myCmpExpvs(const SmallExponent_t* v1, const SmallExponent_t* v2) const
  {
    long deg1=0, deg2=0;
    for (long i=0; i < myNumIndets; ++i)
    {
      deg1 += v1[i];
      deg2 += v2[i];
    }
    if (deg1 != deg2)
    {
      if (deg1 > deg2) return 1; else return -1;
    }

    for (long i = myNumIndets-1; i>0; --i)
      if (v1[i] != v2[i])
      {
        if (v1[i] > v2[i]) return -1; else return 1;
      }
    return 0;
  }


  int PPMonoidEvSmallExpImpl::StdDegRevLexImpl::myCmpWDegPartial(const SmallExponent_t* v1, const SmallExponent_t* v2, long PartialGrDim) const
  {
    CoCoA_ASSERT(PartialGrDim==0 || PartialGrDim==1);
    if (PartialGrDim==0) return 0; // ie: wrt to 0-grading
    long deg1=0, deg2=0;
    for (long i=0; i < myNumIndets; ++i)
    {
      deg1 += v1[i];
      deg2 += v2[i];
    }
    if (deg1 != deg2)
    {
      if (deg1 > deg2) return 1; else return -1;
    }
    return 0;
  }


//--- MatrixOrderingImpl

  PPMonoidEvSmallExpImpl::MatrixOrderingImpl::MatrixOrderingImpl(long NumIndets, long GradingDim, const ConstMatrixView& OrderMatrix):
    CmpBase(NumIndets, GradingDim)
  {
    CoCoA_ASSERT(NumIndets > 0);
    CoCoA_ASSERT(0 <= GradingDim && GradingDim <= NumIndets);
    CoCoA_ASSERT(NumRows(OrderMatrix) == NumIndets);
    CoCoA_ASSERT(NumCols(OrderMatrix) == NumIndets);

    myOrderMatrix.resize(NumIndets, vector<BigInt>(NumIndets));
    for (long i=0; i < NumIndets; ++i)
      for (long j=0; j < NumIndets; ++j)
        myOrderMatrix[i][j] = ConvertTo<BigInt>(OrderMatrix(i,j));
  }


  void PPMonoidEvSmallExpImpl::MatrixOrderingImpl::myWDeg(degree& d, const SmallExponent_t* v) const
  {
    BigInt deg;
    for (long g=0; g < myGradingDim; ++g)
    {
      deg = 0;
      for (long i=0; i < myNumIndets; ++i)
        deg += v[i] * myOrderMatrix[g][i];
      SetComponent(d, g, deg);
    }
  }


  static int MatrixOrderingCmpArrays(const vector<vector<BigInt> >& M, const SmallExponent_t* v1, const SmallExponent_t* v2, long UpTo)
  {
    CoCoA_ASSERT(0 <= UpTo && UpTo <= len(M));
    BigInt s1, s2; // automatically set to 0
    const long ncols = len(M[0]);
    for (long r=0; r < UpTo; ++r)
    {
      for (long i=0; i < ncols; ++i)
      {
        s1 += v1[i]*M[r][i];
        s2 += v2[i]*M[r][i];
      }
      if (s1 != s2)
      {
        if (s1 > s2) return 1; else return -1;
      }
      s1 = 0;
      s2 = 0;
    }
    return 0;
  }


  int PPMonoidEvSmallExpImpl::MatrixOrderingImpl::myCmpExpvs(const SmallExponent_t* v1, const SmallExponent_t* v2) const
  {
    return MatrixOrderingCmpArrays(myOrderMatrix, v1, v2, myNumIndets);
  }


//   int PPMonoidEvSmallExpImpl::MatrixOrderingImpl::myCmpWDeg(const SmallExponent_t* v1, const SmallExponent_t* v2) const
//   {
//     return MatrixOrderingCmpArrays(myOrderMatrix, v1, v2, myGradingDim);
//   }


  int PPMonoidEvSmallExpImpl::MatrixOrderingImpl::myCmpWDegPartial(const SmallExponent_t* v1, const SmallExponent_t* v2, long PartialGrDim) const
  {
    return MatrixOrderingCmpArrays(myOrderMatrix, v1, v2, PartialGrDim);
  }

  //////////////////////////////////////////////////////////////////
  // BIG EXPONENTS

  class PPMonoidEvBigExpImpl: public PPMonoidBase
  {
    typedef PPMonoidElemRawPtr RawPtr;           // just to save typing
    typedef PPMonoidElemConstRawPtr ConstRawPtr; // just to save typing

    class CmpBase
    {
    public:
      CmpBase(long NumIndets, long GradingDim);
      virtual ~CmpBase() {};
    public:
      virtual int myCmpExpvs(const BigInt* v1, const BigInt* v2) const =0;
      virtual void myWDeg(degree& d, const BigInt* v) const =0;
      virtual int myCmpWDegPartial(const BigInt* v1, const BigInt* v2, long GrDim) const =0;
      //      virtual int myCmpWDeg(const BigInt* v1, const BigInt* v2) const =0;
      int myCmpWDeg(const BigInt* v1, const BigInt* v2) const { return myCmpWDegPartial(v1, v2, myGradingDim); }
    protected:
      long myNumIndets;
      long myGradingDim;
    };

    class LexImpl: public CmpBase
    {
    public:
      LexImpl(long NumIndets);
      virtual int myCmpExpvs(const BigInt* v1, const BigInt* v2) const;
      virtual void myWDeg(degree& d, const BigInt* v) const;
      virtual int myCmpWDegPartial(const BigInt* v1, const BigInt* v2, long GrDim) const;
    };


    class StdDegLexImpl: public CmpBase
    {
    public:
      StdDegLexImpl(long NumIndets);
      virtual int myCmpExpvs(const BigInt* v1, const BigInt* v2) const;
      virtual void myWDeg(degree& d, const BigInt* v) const;
      virtual int myCmpWDegPartial(const BigInt* v1, const BigInt* v2, long GrDim) const;
    };


    class StdDegRevLexImpl: public CmpBase
    {
    public:
      StdDegRevLexImpl(long NumIndets);
      virtual int myCmpExpvs(const BigInt* v1, const BigInt* v2) const;
      virtual void myWDeg(degree& d, const BigInt* v) const;
      virtual int myCmpWDegPartial(const BigInt* v1, const BigInt* v2, long GrDim) const;
    };


    class MatrixOrderingImpl: public CmpBase
    {
    public:
      MatrixOrderingImpl(long NumIndets, long GradingDim, const ConstMatrixView& OrderMatrix);
      virtual int myCmpExpvs(const BigInt* v1, const BigInt* v2) const;
      virtual void myWDeg(degree& d, const BigInt* v) const;
      virtual int myCmpWDegPartial(const BigInt* v1, const BigInt* v2, long GrDim) const;
    private:
      std::vector< std::vector<BigInt> > myOrderMatrix;
    };


  public:
    PPMonoidEvBigExpImpl(const std::vector<symbol>& IndetNames, const PPOrdering& ord);
    ~PPMonoidEvBigExpImpl();
  private: // disable copy ctor and assignment
    explicit PPMonoidEvBigExpImpl(const PPMonoidEvBigExpImpl& copy);  // NEVER DEFINED -- copy ctor disabled
    PPMonoidEvBigExpImpl& operator=(const PPMonoidEvBigExpImpl& rhs); // NEVER DEFINED -- assignment disabled

  public:
    void contents() const; // FOR DEBUGGING ONLY

    const std::vector<PPMonoidElem>& myIndets() const;               ///< std::vector whose n-th entry is n-th indet as PPMonoidElem

    // The functions below are operations on power products owned by PPMonoidEvBigExpImpl
    const PPMonoidElem& myOne() const;
    PPMonoidElemRawPtr myNew() const;                                ///< ctor from nothing
    PPMonoidElemRawPtr myNew(PPMonoidElemConstRawPtr rawpp) const;   ///< ctor by assuming ownership
    PPMonoidElemRawPtr myNew(const std::vector<long>& expv) const;   ///< ctor from exp vector
    PPMonoidElemRawPtr myNew(const std::vector<BigInt>& EXPV) const; ///< ctor from exp vector
    void myDelete(RawPtr rawpp) const;                               ///< dtor, frees pp
    void mySwap(RawPtr rawpp1, RawPtr rawpp2) const;                 ///< swap(pp1, pp2);
    void myAssignOne(RawPtr rawpp) const;                            ///< pp = 1
    void myAssign(RawPtr rawpp, ConstRawPtr rawpp1) const;           ///< p = pp1
    void myAssign(RawPtr rawpp, const std::vector<long>& expv) const;///< pp = expv (assign from exp vector)

    void myMul(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const;   ///< pp = pp1*pp2
    using PPMonoidBase::myMulIndetPower;    // disable warnings of overloading
    void myMulIndetPower(RawPtr rawpp, long indet, long exp) const;           ///< pp *= indet^exp, assumes exp >= 0
    void myMulIndetPower(RawPtr rawpp, long indet, const BigInt& EXP) const;  ///< pp *= indet^EXP
    void myDiv(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const;   ///< pp = pp1/pp2
    void myColon(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const; ///< pp = pp1/gcd(pp1,pp2)
    void myGcd(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const;   ///< pp = gcd(pp1,pp2)
    void myLcm(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const;   ///< pp = lcm(pp1,pp2)
    void myRadical(RawPtr rawpp, ConstRawPtr rawpp1) const;                   ///< pp = radical(pp1)
    void myPowerSmallExp(RawPtr rawpp, ConstRawPtr rawpp1, long exp) const;   ///< pp = pp1^exp, assumes exp >= 0
    void myPowerOverflowCheck(ConstRawPtr rawpp1, long exp) const;            ///< throw if pp1^exp would overflow, assumes exp >= 0


    bool myIsOne(ConstRawPtr rawpp) const;                                ///< is pp = 1?
    bool myIsIndet(long& index, ConstRawPtr rawpp) const;                 ///< true iff pp is an indet
    bool myIsIndetPosPower(long& indet, BigInt& EXP, ConstRawPtr rawpp) const;
    bool myIsIndetPosPower(long& indet, long& pow, ConstRawPtr rawpp) const;
    bool myIsIndetPosPower(ConstRawPtr rawpp) const;
    bool myIsCoprime(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const;       ///< are pp1 & pp2 coprime?
    bool myIsEqual(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const;         ///< is pp1 equal to pp2?
    bool myIsDivisible(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const;     ///< does pp2 divide pp1?
    bool myIsSqFree(ConstRawPtr rawpp) const;                             ///< is pp equal to its radical?

    int myCmp(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const;              ///< -1,0,1 as pp1 < = > pp2
    long myStdDeg(ConstRawPtr rawpp) const;                               ///< standard degree of pp
    void myWDeg(degree& d, ConstRawPtr rawpp) const;                      ///< d = grading(pp)
    int myCmpWDeg(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const;          ///< <0, =0, >0 as wdeg(pp1) < = > wdeg(pp2)
    int myCmpWDegPartial(ConstRawPtr rawpp1, ConstRawPtr rawpp2, long GrDim ) const; ///< as myCmpWDeg wrt the first GrDim weights
    long myExponent(ConstRawPtr rawpp, long indet) const;                 ///< exponent of indet in pp
    void myBigExponent(BigInt& EXP, ConstRawPtr rawpp, long indet) const; ///< EXP = exponent of indet in pp
    void myExponents(std::vector<long>& expv, ConstRawPtr rawpp) const;   ///< expv[i] = exponent(pp,i)
    void myBigExponents(std::vector<BigInt>& v, ConstRawPtr rawpp) const; ///< get exponents, SHOULD BE myExponents ???
    void myIndetsIn(std::vector<bool>& v, ConstRawPtr rawpp) const;       ///< v[i] = true if i-th indet has exponent != 0
    void myOutputSelf(std::ostream& out) const;                           ///< out << PPM
    // INHERITED DEFINITION of virtual  void myOutput(std::ostream& out, ConstRawPtr rawpp) const;
    void myDebugPrint(std::ostream& out, ConstRawPtr rawpp) const;        ///< print pp in debugging format???


  private: // auxiliary functions
    BigInt* myExpv(RawPtr) const;
    const BigInt* myExpv(ConstRawPtr) const;

    void myComputeDivMask(DivMask& dm, const DivMaskRule& DivMaskImpl, ConstRawPtr rawpp) const; ///< used by PPWithMask
    bool myCheckExponents(const std::vector<long>& expv) const;
    bool myCheckExponents(const std::vector<BigInt>& EXPV) const;

  private: // data members
    ///@name Class members
    //@{
    const long myEntrySize;     ///< size in ????
    mutable MemPool myMemMgr;     // IMPORTANT: this must come *before* myIndetVector and myOnePtr.
    std::vector<PPMonoidElem> myIndetVector; ///< the indets as PPMonoidElems
    std::unique_ptr<CmpBase> myOrdPtr;         ///< actual implementation of the ordering [should be const???]
    std::unique_ptr<PPMonoidElem> myOnePtr;
    //@}
  };


  // File local inline functions

  inline BigInt* PPMonoidEvBigExpImpl::myExpv(RawPtr rawpp) const
  {
    return static_cast<BigInt*>(rawpp.myRawPtr());
  }


  inline const BigInt* PPMonoidEvBigExpImpl::myExpv(ConstRawPtr rawpp) const
  {
    return static_cast<const BigInt*>(rawpp.myRawPtr());
  }


  bool PPMonoidEvBigExpImpl::myCheckExponents(const std::vector<long>& expv) const
  {
    // Check len(expv) == myNumIndets.
    // Check exps are non-neg and not too big.
    if (len(expv) != myNumIndets) return false;
    for (long i=0; i < myNumIndets; ++i)
      if (expv[i] < 0) return false;
    return true;
  }

  bool PPMonoidEvBigExpImpl::myCheckExponents(const std::vector<BigInt>& expv) const
  {
    // Check len(expv) == myNumIndets.
    // Check exps are non-neg and not too big.
    if (len(expv) != myNumIndets) return false;
    for (long i=0; i < myNumIndets; ++i)
      if (expv[i] < 0) return false;
    return true;
  }


  //----   Constructors & destructor   ----//

  PPMonoidEvBigExpImpl::PPMonoidEvBigExpImpl(const std::vector<symbol>& IndetNames, const PPOrdering& ord):
      PPMonoidBase(ord, IndetNames),
      myEntrySize(sizeof(BigInt)*NumIndets(ord)),
      myMemMgr(myEntrySize, "PPMonoidEvBigExpImpl.myMemMgr"),
      myIndetVector()
  {
    // Put the correct implementation in myOrdPtr...  [VERY UGLY CODE!!!]
    if (IsLex(ord))
      myOrdPtr.reset(new LexImpl(myNumIndets));
    else if (IsStdDegLex(ord))
      myOrdPtr.reset(new StdDegLexImpl(myNumIndets));
    else if (IsStdDegRevLex(ord))
      myOrdPtr.reset(new StdDegRevLexImpl(myNumIndets));
    else
      myOrdPtr.reset(new MatrixOrderingImpl(myNumIndets, GradingDim(ord), OrdMat(ord)));

    myRefCountInc();  // this is needed for exception cleanliness, in case one of the lines below throws
    myOnePtr.reset(new PPMonoidElem(PPMonoid(this)));
    {
      // IMPORTANT: this block destroys pp *before* the call to myRefCountZero.
      PPMonoidElem pp(PPMonoid(this));
      vector<long> expv(myNumIndets);
      for (long i=0; i < myNumIndets; ++i)
      {
        expv[i] = 1;
        myAssign(raw(pp), expv);
        myIndetVector.push_back(pp);
        expv[i] = 0;
      }
    }
    myRefCountZero();
  }


  PPMonoidEvBigExpImpl::~PPMonoidEvBigExpImpl()
  {}

/////////////////////////////////////////////////////////////////////////////


  const std::vector<PPMonoidElem>& PPMonoidEvBigExpImpl::myIndets() const
  {
    return myIndetVector;
  }


  const PPMonoidElem& PPMonoidEvBigExpImpl::myOne() const
  {
    return *myOnePtr;
  }


  PPMonoidElemRawPtr PPMonoidEvBigExpImpl::myNew() const
  {
    PPMonoidElemRawPtr rawpp(myMemMgr.alloc());
    BigInt* const expv = myExpv(rawpp);
    for (long i=0; i < myNumIndets; ++i) // BUG BUG: if this throws, rawpp is leaked
      new(&expv[i]) BigInt;
    myAssignOne(rawpp); // cannot throw
    return rawpp;
  }

  PPMonoidElemRawPtr PPMonoidEvBigExpImpl::myNew(PPMonoidElemConstRawPtr rawcopypp) const
  {
    PPMonoidElemRawPtr rawpp(myMemMgr.alloc());
    BigInt* const expv = myExpv(rawpp);
    for (long i=0; i < myNumIndets; ++i) // BUG BUG: if this throws, rawpp is leaked
      new(&expv[i]) BigInt;
    myAssign(rawpp, rawcopypp); // cannot throw
    return rawpp;
  }


  PPMonoidElemRawPtr PPMonoidEvBigExpImpl::myNew(const std::vector<long>& expv) const
  {
    CoCoA_ASSERT(myCheckExponents(expv));
    PPMonoidElemRawPtr rawpp(myMemMgr.alloc());
    BigInt* const expv1 = myExpv(rawpp);
    for (long i=0; i < myNumIndets; ++i) // BUG BUG: if this throws, rawpp is leaked
      new(&expv1[i]) BigInt;
    myAssign(rawpp, expv); // cannot throw
    return rawpp;
  }


  PPMonoidElemRawPtr PPMonoidEvBigExpImpl::myNew(const std::vector<BigInt>& EXPV) const
  {
    CoCoA_ASSERT(myCheckExponents(EXPV));
    PPMonoidElemRawPtr rawpp(myMemMgr.alloc());
    BigInt* const expv = myExpv(rawpp);
    for (long i=0; i < myNumIndets; ++i) // BUG BUG: if this throws, rawpp is leaked
      new(&expv[i]) BigInt;
    for (long i = 0; i < myNumIndets; ++i)
      expv[i] = EXPV[i];
    return rawpp;
  }


  void PPMonoidEvBigExpImpl::myAssignOne(RawPtr rawpp) const
  {
    BigInt* const expv = myExpv(rawpp);
    for (long i = 0; i < myNumIndets; ++i)
      expv[i] = 0;
  }


  void PPMonoidEvBigExpImpl::myAssign(RawPtr rawpp, ConstRawPtr rawpp1) const
  {
    if (rawpp == rawpp1) return;

    BigInt* const expv = myExpv(rawpp);
    const BigInt* const expv1 = myExpv(rawpp1);
    std::copy(&expv1[0], &expv1[myNumIndets], &expv[0]);
  }

  void PPMonoidEvBigExpImpl::myAssign(RawPtr rawpp, const vector<long>& expv1) const
  {
    CoCoA_ASSERT(myCheckExponents(expv1));

    BigInt* const expv = myExpv(rawpp);
    for (long i = 0; i < myNumIndets; ++i)
      expv[i] = expv1[i];
  }


  void PPMonoidEvBigExpImpl::myDelete(RawPtr rawpp) const
  {
    BigInt* const expv = myExpv(rawpp);;
    for (long i=0; i < myNumIndets; ++i)
      expv[i].~BigInt();
    myMemMgr.free(rawpp.myRawPtr());
  }


  void PPMonoidEvBigExpImpl::mySwap(RawPtr rawpp1, RawPtr rawpp2) const
  {
    if (rawpp1 == rawpp2) return;
    BigInt* const expv1 = myExpv(rawpp1);
    BigInt* expv2 = myExpv(rawpp2);
    for (long i = 0; i < myNumIndets; ++i)
      std::swap(expv1[i], expv2[i]);
  }


  void PPMonoidEvBigExpImpl::myMul(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    BigInt* const expv = myExpv(rawpp);
    const BigInt* const expv1 = myExpv(rawpp1);
    const BigInt* const expv2 = myExpv(rawpp2);
    for (long i = 0; i < myNumIndets; ++i)
      expv[i] = expv1[i] + expv2[i];
  }


  void PPMonoidEvBigExpImpl::myMulIndetPower(RawPtr rawpp, long indet, long exp) const  // assumes exp >= 0
  {
    CoCoA_ASSERT(exp >= 0);
    CoCoA_ASSERT(0 <= indet && indet < myNumIndets);
    BigInt* const expv = myExpv(rawpp);
    expv[indet] += exp;
  }


  // This fn overrides the default defn in PPMonoid.C (which throws for big exps)
  void PPMonoidEvBigExpImpl::myMulIndetPower(RawPtr rawpp, long indet, const BigInt& EXP) const
  {
    CoCoA_ASSERT(0 <= indet && indet < myNumIndets);
    CoCoA_ASSERT(EXP >= 0);
    BigInt* const expv = myExpv(rawpp);
    expv[indet] += EXP;
  }


  void PPMonoidEvBigExpImpl::myDiv(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    BigInt* const expv = myExpv(rawpp);
    const BigInt* const expv1 = myExpv(rawpp1);
    const BigInt* const expv2 = myExpv(rawpp2);
    for (long i = 0; i < myNumIndets; ++i)
      expv[i] = expv1[i] - expv2[i]; // assert that result is positive?
  }


  void PPMonoidEvBigExpImpl::myColon(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    // No worries about aliasing.
    BigInt* const expv = myExpv(rawpp);
    const BigInt* const expv1 = myExpv(rawpp1);
    const BigInt* const expv2 = myExpv(rawpp2);

    for (long i = 0; i < myNumIndets; ++i)
      if (expv1[i] > expv2[i])
        expv[i] = expv1[i] - expv2[i];
      else
        expv[i] = 0;
  }


  void PPMonoidEvBigExpImpl::myGcd(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    // No worries about aliasing.
    BigInt* const expv = myExpv(rawpp);
    const BigInt* const expv1 = myExpv(rawpp1);
    const BigInt* const expv2 = myExpv(rawpp2);

    for (long i = 0; i < myNumIndets; ++i)
      expv[i] = min(expv1[i], expv2[i]);
  }


  void PPMonoidEvBigExpImpl::myLcm(RawPtr rawpp, ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    // No worries about aliasing.
    BigInt* const expv = myExpv(rawpp);
    const BigInt* const expv1 = myExpv(rawpp1);
    const BigInt* const expv2 = myExpv(rawpp2);

    for (long i = 0; i < myNumIndets; ++i)
      expv[i] = max(expv1[i], expv2[i]);
  }


  void PPMonoidEvBigExpImpl::myRadical(RawPtr rawpp, ConstRawPtr rawpp1) const
  {
    BigInt* const expv = myExpv(rawpp);
    const BigInt* const expv1 = myExpv(rawpp1);

    for (long i = 0; i < myNumIndets; ++i)
      expv[i] = (expv1[i] > 0);
  }


  void PPMonoidEvBigExpImpl::myPowerSmallExp(RawPtr rawpp, ConstRawPtr rawpp1, long exp) const  // assumes exp >= 0
  {
    CoCoA_ASSERT(exp >= 0);
    BigInt* const expv = myExpv(rawpp);
    const BigInt* const expv1 = myExpv(rawpp1);

    for (long i = 0; i < myNumIndets; ++i)
      expv[i] = exp * expv1[i];
  }


  void PPMonoidEvBigExpImpl::myPowerOverflowCheck(ConstRawPtr /*rawpp*/, long exp) const
  {
    (void)(exp); // just to avoid compiler warning
    CoCoA_ASSERT(exp >= 0);
    // nothing to check -- exps here cannot overflow :-)
  }


  bool PPMonoidEvBigExpImpl::myIsOne(ConstRawPtr rawpp) const
  {
    const BigInt* const expv = myExpv(rawpp);

    for (long i = 0; i < myNumIndets; ++i)
      if (!IsZero(expv[i])) return false;

    return true;
  }


  bool PPMonoidEvBigExpImpl::myIsIndet(long& index, ConstRawPtr rawpp) const
  {
    const BigInt* const expv = myExpv(rawpp);
    long j = myNumIndets;
    for (long i = 0; i < myNumIndets; ++i)
    {
      if (IsZero(expv[i])) continue;
      if (j != myNumIndets || !IsOne(expv[i])) return false;
      j = i;
    }
    if (j == myNumIndets) return false;
    index = j;
    return true;
  }


  bool PPMonoidEvBigExpImpl::myIsIndetPosPower(long& index, BigInt& EXP, ConstRawPtr rawpp) const
  {
    const BigInt* const expv = myExpv(rawpp);
    long j = myNumIndets;
    for (long i = 0; i < myNumIndets; ++i)
    {
      if (IsZero(expv[i])) continue;
      if (j != myNumIndets) return false;
      j = i;
    }
    if (j == myNumIndets) return false;
    index = j;
    EXP = expv[index];
    return true;
  }

  bool PPMonoidEvBigExpImpl::myIsIndetPosPower(long& index, long& pow, ConstRawPtr rawpp) const
  {
    BigInt POW;
    long TmpIndex, TmpPow;
    if (!myIsIndetPosPower(TmpIndex, POW, rawpp)) return false;
    if (!IsConvertible(TmpPow, POW))
      CoCoA_ERROR(ERR::ArgTooBig, "PPMonoidEvBigExpImpl::myIsIndetPosPower");
    index = TmpIndex;
    pow = TmpPow;
    return true;
  }

  bool PPMonoidEvBigExpImpl::myIsIndetPosPower(ConstRawPtr rawpp) const
  {
    const BigInt* const expv = myExpv(rawpp);
    long j = myNumIndets;
    for (long i = 0; i < myNumIndets; ++i)
    {
      if (IsZero(expv[i])) continue;
      if (j != myNumIndets) return false;
      j = i;
    }
    return j != myNumIndets;
  }


  bool PPMonoidEvBigExpImpl::myIsCoprime(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    const BigInt* const expv1 = myExpv(rawpp1);
    const BigInt* const expv2 = myExpv(rawpp2);

    for (long i = 0; i < myNumIndets; ++i)
      if (!IsZero(expv1[i]) && !IsZero(expv2[i])) return false;

    return true;
  }


  bool PPMonoidEvBigExpImpl::myIsEqual(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    const BigInt* const expv1 = myExpv(rawpp1);
    const BigInt* const expv2 = myExpv(rawpp2);

    for (long i = 0; i < myNumIndets; ++i)
      if (expv1[i] != expv2[i]) return false;
    return true;
  }


  bool PPMonoidEvBigExpImpl::myIsDivisible(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    const BigInt* const expv1 = myExpv(rawpp1);
    const BigInt* const expv2 = myExpv(rawpp2);

    for (long i = 0; i < myNumIndets; ++i)
      if (expv1[i] < expv2[i]) return false;

    return true;
  }


  bool PPMonoidEvBigExpImpl::myIsSqFree(ConstRawPtr rawpp) const
  {
    const BigInt* const expv = myExpv(rawpp);

    for (long i = 0; i < myNumIndets; ++i)
      if (expv[i] > 1) return false;

    return true;
  }


  int PPMonoidEvBigExpImpl::myCmp(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    return myOrdPtr->myCmpExpvs(myExpv(rawpp1), myExpv(rawpp2));
  }


  long PPMonoidEvBigExpImpl::myStdDeg(ConstRawPtr rawpp) const
  {
    const BigInt* const expv = myExpv(rawpp);
    BigInt d;
    for (long i=0; i < myNumIndets; ++i)
      d += expv[i];
    long ans;
    if (!IsConvertible(ans, d))
      CoCoA_ERROR(ERR::ArgTooBig, "PPMonoidEvBigExpImpl::myStdDeg");
    return ans;
  }


  void PPMonoidEvBigExpImpl::myWDeg(degree& d, ConstRawPtr rawpp) const
  {
    myOrdPtr->myWDeg(d, myExpv(rawpp));
  }


  int PPMonoidEvBigExpImpl::myCmpWDeg(ConstRawPtr rawpp1, ConstRawPtr rawpp2) const
  {
    return myOrdPtr->myCmpWDeg(myExpv(rawpp1), myExpv(rawpp2));
  }


  int PPMonoidEvBigExpImpl::myCmpWDegPartial(ConstRawPtr rawpp1, ConstRawPtr rawpp2, long GrDim) const
  {
    return myOrdPtr->myCmpWDegPartial(myExpv(rawpp1), myExpv(rawpp2), GrDim);
  }


  long PPMonoidEvBigExpImpl::myExponent(ConstRawPtr rawpp, long indet) const
  {
    CoCoA_ASSERT(0 <= indet && indet < myNumIndets);
    long ans;
    if (!IsConvertible(ans, myExpv(rawpp)[indet]))
      CoCoA_ERROR(ERR::ArgTooBig, "PPMonoidEvBigExpImpl::myExponent");
    return ans;
  }

  void PPMonoidEvBigExpImpl::myBigExponent(BigInt& EXP, ConstRawPtr rawpp, long indet) const
  {
    CoCoA_ASSERT(0 <= indet && indet < myNumIndets);
    EXP = myExpv(rawpp)[indet];
  }


  void PPMonoidEvBigExpImpl::myExponents(vector<long>& v, ConstRawPtr rawpp) const
  {
    CoCoA_ASSERT(len(v) == myNumIndets);
    const BigInt* const expv = myExpv(rawpp);
    bool OK = true;
    for (long i=0; i < myNumIndets; ++i)
    {
      if (!IsConvertible(v[i], expv[i]))
      {
        OK = false;
        if (expv[i] > 0)
          v[i] = numeric_limits<long>::max();
        else
          v[i] = numeric_limits<long>::min();
      }
      if (!OK)
        CoCoA_ERROR(ERR::ArgTooBig, "PPMonoidEvBigExpImpl::myExponents");
    }
  }


  void PPMonoidEvBigExpImpl::myBigExponents(vector<BigInt>& expv, ConstRawPtr rawpp) const
  {
    CoCoA_ASSERT(len(expv) == myNumIndets);
    const BigInt* const v = myExpv(rawpp);
    for (long i=0; i < myNumIndets; ++i)  expv[i] = v[i];
  }


  void PPMonoidEvBigExpImpl::myIndetsIn(std::vector<bool>& v, ConstRawPtr rawpp) const
  {
    CoCoA_ASSERT(len(v) == myNumIndets);
    const BigInt* const expv = myExpv(rawpp);
    for (long i=0; i < myNumIndets; ++i)
      if (!IsZero(expv[i]))
        v[i] = true;
  }


  void PPMonoidEvBigExpImpl::myComputeDivMask(DivMask& /*dm*/, const DivMaskRule& /*DivMaskImpl*/, ConstRawPtr /*rawpp*/) const
  {
    CoCoA_ERROR(ERR::NYI, "PPMonoidEvBigExpImpl::myComputeDivMask");
// BUG BUG can't do this (yet)    DivMaskImpl->myAssignFromExpv(dm, myExpv(rawpp), myNumIndets);
  }


  void PPMonoidEvBigExpImpl::myOutputSelf(std::ostream& out) const
  {
    if (!out) return;  // short-cut for bad ostreams
    out << "PPMonoidBigEv(" << myNumIndets << ", " << myOrd <<")";
  }


  void PPMonoidEvBigExpImpl::myDebugPrint(std::ostream& out, ConstRawPtr rawpp) const
  {
    if (!out) return;  // short-cut for bad ostreams
    out << "DEBUG PP: myNumIndets=" << myNumIndets << ", exps=[";
    for (long i=0; i < myNumIndets; ++i)
      out << myExponent(rawpp, i) << " ";
    out << "]" << std::endl;
  }


//--- orderings ----------------------------------------

  PPMonoidEvBigExpImpl::CmpBase::CmpBase(long NumIndets, long GradingDim):
    myNumIndets(NumIndets),
    myGradingDim(GradingDim)
  {
    CoCoA_ASSERT(NumIndets > 0);
    CoCoA_ASSERT(NumIndets < 1000000); // complain about ridiculously large number of indets
    CoCoA_ASSERT(0 <= GradingDim && GradingDim <= NumIndets);
  }


//--- LexImpl

  PPMonoidEvBigExpImpl::LexImpl::LexImpl(long NumIndets):
    CmpBase(NumIndets, 0)
  {
  }


  int PPMonoidEvBigExpImpl::LexImpl::myCmpExpvs(const BigInt* v1, const BigInt* v2) const
  {
    for (long i=0; i<myNumIndets; ++i)
      if (v1[i] != v2[i]) return (v1[i]>v2[i] ? 1 : -1 );
    return 0;
  }


  void PPMonoidEvBigExpImpl::LexImpl::myWDeg(degree& /*d*/, const BigInt* /*v*/) const
  {
    // deliberately does nothing because GradingDim=0
    //??? should it assign:  d = []
  }


  int PPMonoidEvBigExpImpl::LexImpl::myCmpWDegPartial(const BigInt* /*v1*/, const BigInt* /*v2*/, long /*GrDim*/) const
  {
    return 0; // GradingDim=0, and all degrees in Z^0 are equal
  }

//--- StdDegLexImpl

  PPMonoidEvBigExpImpl::StdDegLexImpl::StdDegLexImpl(long NumIndets):
    CmpBase(NumIndets, 1)
  {
  }


  void PPMonoidEvBigExpImpl::StdDegLexImpl::myWDeg(degree& d, const BigInt* v) const
  {
    BigInt deg;
    for (long i=0; i < myNumIndets; ++i)
      deg += v[i];
    SetComponent(d, 0, deg);
  }


  int PPMonoidEvBigExpImpl::StdDegLexImpl::myCmpExpvs(const BigInt* v1, const BigInt* v2) const
  {
    BigInt deg1;
    BigInt deg2;
    for (long i=0; i < myNumIndets; ++i)
    {
      deg1 += v1[i];
      deg2 += v2[i];
    }
    if (deg1 != deg2) return (deg1>deg2 ? 1 : -1 );

    for (long i=0; i < myNumIndets; ++i)
      if (v1[i] != v2[i]) return (v1[i]>v2[i] ? 1 : -1 );
    return 0;
  }


  int PPMonoidEvBigExpImpl::StdDegLexImpl::myCmpWDegPartial(const BigInt* v1, const BigInt* v2, long GrDim) const
  {
    CoCoA_ASSERT(0 <= GrDim && GrDim <= 1);
    if (GrDim == 0) return 0; // ie: wrt to 0-grading
    BigInt deg1;
    BigInt deg2;
    for (long i=0; i < myNumIndets; ++i)
    {
      deg1 += v1[i];
      deg2 += v2[i];
    }
    if (deg1 != deg2) return (deg1>deg2 ? 1 : -1 );
    return 0;
  }


//--- StdDegRevLexImpl

  PPMonoidEvBigExpImpl::StdDegRevLexImpl::StdDegRevLexImpl(long NumIndets):
    CmpBase(NumIndets, 1)
  {
  }


  void PPMonoidEvBigExpImpl::StdDegRevLexImpl::myWDeg(degree& d, const BigInt* v) const
  {
    BigInt deg;
    for (long i=0; i < myNumIndets; ++i)
      deg += v[i];
    SetComponent(d, 0, deg);
  }


  int PPMonoidEvBigExpImpl::StdDegRevLexImpl::myCmpExpvs(const BigInt* v1, const BigInt* v2) const
  {
    BigInt deg1;
    BigInt deg2;
    for (long i=0; i < myNumIndets; ++i)
    {
      deg1 += v1[i];
      deg2 += v2[i];
    }
    if (deg1 != deg2) return (deg1>deg2 ? 1 : -1 );

    for (long i = myNumIndets-1; i>0; --i)
      if (v1[i] != v2[i]) return (v1[i]>v2[i] ? -1 : 1 );
    return 0;
  }


  int PPMonoidEvBigExpImpl::StdDegRevLexImpl::myCmpWDegPartial(const BigInt* v1, const BigInt* v2, long GrDim) const
  {
    CoCoA_ASSERT(0 <= GrDim && GrDim <= 1);
    if (GrDim == 0) return 0; // ie: wrt to 0-grading
    BigInt deg1;
    BigInt deg2;
    for (long i=0; i < myNumIndets; ++i)
    {
      deg1 += v1[i];
      deg2 += v2[i];
    }
    if (deg1 != deg2) return (deg1>deg2 ? 1 : -1 );
    return 0;
  }


//--- MatrixOrderingImpl

  PPMonoidEvBigExpImpl::MatrixOrderingImpl::MatrixOrderingImpl(long NumIndets, long GradingDim, const ConstMatrixView& OrderMatrix):
    CmpBase(NumIndets, GradingDim)
  {
    CoCoA_ASSERT(0 <= NumIndets);
    CoCoA_ASSERT(0 <= GradingDim && GradingDim <= NumIndets);
    CoCoA_ASSERT(NumRows(OrderMatrix) == NumIndets);
    CoCoA_ASSERT(NumCols(OrderMatrix) == NumIndets);

    myOrderMatrix.resize(NumIndets, vector<BigInt>(NumIndets));
    for (long i=0; i < NumIndets; ++i)
      for (long j=0; j < NumIndets; ++j)
        myOrderMatrix[i][j] = ConvertTo<BigInt>(OrderMatrix(i,j));
  }


  void PPMonoidEvBigExpImpl::MatrixOrderingImpl::myWDeg(degree& d, const BigInt* v) const
  {
    BigInt deg;
    for (long g=0; g < myGradingDim; ++g)
    {
      deg = 0;
      for (long i=0; i < myNumIndets; ++i)
        deg += v[i] * myOrderMatrix[g][i];
      SetComponent(d, g, deg);
    }
  }


  static int MatrixOrderingCmpArrays(const vector<vector<BigInt> >& M, const BigInt* v1, const BigInt* v2, long UpTo)
  {
    CoCoA_ASSERT(0 <= UpTo && UpTo <= len(M));
    BigInt s1, s2; // automatically set to 0
    const long ncols = len(M[0]);
    for (long r=0; r < UpTo; ++r)
    {
      for (long i=0; i < ncols; ++i)
      {
        s1 += v1[i]*M[r][i];
        s2 += v2[i]*M[r][i];
      }
      if (s1 != s2) return (s1>s2 ? 1 : -1 );
      s1 = 0;
      s2 = 0;
    }
    return 0;
  }


  int PPMonoidEvBigExpImpl::MatrixOrderingImpl::myCmpExpvs(const BigInt* v1, const BigInt* v2) const
  {
    return MatrixOrderingCmpArrays(myOrderMatrix, v1, v2, myNumIndets);
  }


  int PPMonoidEvBigExpImpl::MatrixOrderingImpl::myCmpWDegPartial(const BigInt* v1, const BigInt* v2, long GrDim) const
  {
    CoCoA_ASSERT(0 <= GrDim && GrDim <= myGradingDim);
    return MatrixOrderingCmpArrays(myOrderMatrix, v1, v2, GrDim);
  }


  //////////////////////////////////////////////////////////////////
  // Pseudo-ctors

  PPMonoid NewPPMonoidEv(const std::vector<symbol>& IndetNames, const PPOrdering& ord, PPExpSize ExpSize)
  {
    // Sanity check on the indet names given.
    const long nvars = NumIndets(ord);

    if (len(IndetNames) != nvars)
      CoCoA_ERROR(ERR::BadNumIndets, "NewPPMonoidEv(IndetNames,ord)");
    if (!AreDistinct(IndetNames))
      CoCoA_ERROR(ERR::BadIndetNames, "NewPPMonoidEv(IndetNames,ord)");
    if (!AreArityConsistent(IndetNames))
      CoCoA_ERROR(ERR::BadIndetNames, "NewPPMonoidEv(IndetNames,ord)");

    if (ExpSize == SmallExps)
      return PPMonoid(new PPMonoidEvSmallExpImpl(IndetNames, ord));
    return PPMonoid(new PPMonoidEvBigExpImpl(IndetNames, ord));
  }

  PPMonoid NewPPMonoidEv(const std::vector<symbol>& IndetNames, const PPOrderingCtor& OrdCtor, PPExpSize ExpSize)
  {
    return NewPPMonoidEv(IndetNames, OrdCtor(len(IndetNames)), ExpSize);
  }


} // end of namespace CoCoA


// RCS header/log in the next few lines
// $Header: /Volumes/Home_1/cocoa/cvs-repository/CoCoALib-0.99/src/AlgebraicCore/PPMonoidEv.C,v 1.50 2020/02/15 11:10:22 abbott Exp $
// $Log: PPMonoidEv.C,v $
// Revision 1.50  2020/02/15 11:10:22  abbott
// Summary: IsIndetPosPower now gives false instead of error for input 1
//
// Revision 1.49  2020/02/11 16:56:41  abbott
// Summary: Corrected last update (see redmine 969)
//
// Revision 1.48  2020/02/11 16:12:18  abbott
// Summary: Added some checks for bad ostream (even to mem fns myOutput); see redmine 969
//
// Revision 1.47  2019/03/04 10:29:36  abbott
// Summary: Changed auto_ptr into unqiue_ptr
//
// Revision 1.46  2018/05/18 12:15:04  bigatti
// -- renamed IntOperations --> BigIntOps
//
// Revision 1.45  2017/04/18 12:50:06  abbott
// Summary: Corrected ifdef use of CoCoA_THREADSAFE_HACK and CoCoA_DEBUG
//
// Revision 1.44  2016/11/03 12:25:25  abbott
// Summary: Changed IsRadical (for PPMonoidElem) into IsSqFree
//
// Revision 1.43  2015/12/01 13:11:01  abbott
// Summary: Changed mem fn PPOrderingCtor::myCtor into operator(); also for ModuleOrderingCtor; see issue 829
//
// Revision 1.42  2015/11/30 21:53:55  abbott
// Summary: Major update to matrices for orderings (not yet complete, some tests fail)
//
// Revision 1.41  2015/06/30 12:54:07  abbott
// Summary: Added new fn myIndetsIn
// Author: JAA
//
// Revision 1.40  2015/04/16 20:19:39  abbott
// Summary: Silly change to avoid a compiler warning about unused param
// Author: JAA
//
// Revision 1.39  2015/04/16 16:36:33  abbott
// Summary: Cleaned impls of myPowerOverflowCheck
// Author: JAA
//
// Revision 1.38  2015/04/13 16:17:58  abbott
// Summary: Corrected silly typo.
// Author: JAA
//
// Revision 1.37  2015/04/13 14:42:07  abbott
// Summary: Added myPowerOverflowCheck (1st version)
// Author: JAA
//
// Revision 1.36  2014/07/31 13:10:45  bigatti
// -- GetMatrix(PPO) --> OrdMat(PPO)
// -- added OrdMat and GradingMat to PPOrdering, PPMonoid, SparsePolyRing
//
// Revision 1.35  2014/07/03 15:36:35  abbott
// Summary: Cleaned up impl of PPMonoids: moved myIndetSymbols & myNumIndets to base class
// Author: JAA
//
// Revision 1.34  2014/05/14 15:57:15  bigatti
// -- added "using" for clang with superpedantic flag
//
// Revision 1.33  2014/01/28 09:59:10  abbott
// Replaced two calls to IsInteger by calls to ConvertTo<BigInt>.
//
// Revision 1.32  2014/01/16 16:11:43  abbott
// Added some missing const keywords (to local variables)
//
// Revision 1.31  2013/03/15 11:00:50  abbott
// Added check for exponent overflow when powering a PP.
// Merged PPMonoidEv and PPMonoidEvZZ implementations into a single file.
// Implemented new interface for pseudo-ctors for PPMonoidEv which uses a "flag"
// to say whether exponents are big or not.
//
// Revision 1.30  2012/05/28 09:18:21  abbott
// Created IntOperations which gathers together all operations on
// integers (both big and small).  Many consequential changes.
//
// Revision 1.29  2012/02/10 17:04:23  abbott
// Removed a commented out decl of a member fn: myNew(vector<BigInt>)
//
// Revision 1.28  2012/01/26 16:50:55  bigatti
// -- changed back_inserter into insert
//
// Revision 1.27  2011/08/14 15:52:17  abbott
// Changed ZZ into BigInt (phase 1: just the library sources).
//
// Revision 1.26  2011/05/03 12:13:12  abbott
// Added static const data member ourMaxExp.
// Code is more readable, and compiler doesn't grumble any more.
//
// Revision 1.25  2011/03/22 23:12:16  abbott
// Removed some spurious junk at the start of the file.
// Corrected copyright years.
//
// Revision 1.24  2011/03/22 22:38:15  abbott
// Fixed some wrong static_casts inside some CoCoA_ASSERTs.
//
// Revision 1.23  2011/03/10 17:47:42  bigatti
// -- fixed a CoCoA_ASSERT
//
// Revision 1.22  2011/03/10 17:26:19  bigatti
// -- changed unsigned long into long in soe CoCoA_ASSERT
//
// Revision 1.21  2011/03/10 16:39:34  abbott
// Replaced (very many) size_t by long in function interfaces (for rings,
// PPMonoids and modules).  Also replaced most size_t inside fn defns.
//
// Revision 1.20  2011/03/01 14:14:54  bigatti
// -- made CmpBase dtor non-inline for problems with g++ on some ppc
//
// Revision 1.19  2010/12/17 16:09:51  abbott
// Corrected used of myIndetSymbols in some assertions.
//
// Revision 1.18  2010/11/30 11:18:11  bigatti
// -- renamed IndetName --> IndetSymbol
//
// Revision 1.17  2010/11/05 16:21:08  bigatti
// -- added ZZExponents
//
// Revision 1.16  2010/10/06 14:10:24  abbott
// Added increments to the ref count in ring and PPMonoid ctors to make
// them exception safe.
//
// Revision 1.15  2010/02/03 16:13:52  abbott
// Added new single word tags for specifying the ordering in PPMonoid
// pseudo-ctors.
//
// Revision 1.14  2010/02/02 16:44:31  abbott
// Added radical & IsRadical (via mem fns myRadical & myIsRadical)
// for PPMonoidElems.
//
// Revision 1.13  2009/12/23 18:53:52  abbott
// Major change to conversion functions:
//   convert(..) is now a procedure instead of a function
//   IsConvertible(..) replaces the former convert(..) function
//   Added new NumericCast conversion function (placeholder for BOOST feature)
//   Consequent changes in code which uses these features.
//
// Revision 1.12  2009/09/22 14:01:33  bigatti
// -- added myCmpWDegPartial (ugly name, I know....)
// -- cleaned up and realigned code in PPMonoid*.C files
//
// Revision 1.11  2009/07/02 16:34:18  abbott
// Minor change to keep the compiler quiet.
//
// Revision 1.10  2008/03/26 16:52:04  abbott
// Added exponent overflow checks (also for ordvs) when CoCoA_DEBUG is active.
//
// Revision 1.9  2007/12/05 11:06:24  bigatti
// -- changed "size_t StdDeg/myStdDeg(f)" into "long"  (and related functions)
// -- changed "log/myLog(f, i)" into "MaxExponent/myMaxExponent(f, i)"
// -- fixed bug in "IsOne(ideal)" in SparsePolyRing.C
//
// Revision 1.8  2007/12/04 14:27:07  bigatti
// -- changed "log(pp, i)" into "exponent(pp, i)"
//
// Revision 1.7  2007/10/30 17:14:07  abbott
// Changed licence from GPL-2 only to GPL-3 or later.
// New version for such an important change.
//
// Revision 1.6  2007/09/25 16:32:30  abbott
// Several minor changes to silence gcc-4.3:
//    more #includes,
//    and fixed a template problemm in RegisterServerOps.C
//
// Revision 1.5  2007/05/31 14:54:31  bigatti
// -- now using AreDistinct and AreArityConsistent for sanity check on
//    indet names
//
// Revision 1.3  2007/05/03 10:35:23  abbott
// Added new PPMonoidEvZZ with (virtually) unlimited exponents.
// Modified test-PPMonoid1.C accordingly.
// Added warning in doc about silent/unchecked exponent overflow in other
// PPMonoids.
//
// Revision 1.2  2007/03/23 18:38:42  abbott
// Separated the "convert" functions (and CheckedCast) into their own files.
// Many consequential changes.  Also corrected conversion to doubles.
//
// Revision 1.1.1.1  2007/03/09 15:16:11  abbott
// Imported files
//
// Revision 1.15  2007/03/08 18:22:29  cocoa
// Just whitespace cleaning.
//
// Revision 1.14  2007/03/08 17:43:11  cocoa
// Swapped order of args to the NewPPMonoid pseudo ctors.
//
// Revision 1.13  2007/03/08 11:07:12  cocoa
// Made pseudo ctors for polynomial rings more uniform.  This allowed me to
// remove an include of CoCoA/symbol.H  from the RingDistrM*.H files, but then
// I had to put the include in several .C files.
//
// Revision 1.12  2006/12/07 17:23:46  cocoa
// -- for compilation with _Wextra: commented out names of unused arguments
//
// Revision 1.11  2006/12/06 17:35:58  cocoa
// -- style: RawPtr args are now called "raw.."
//
// Revision 1.10  2006/11/27 13:06:23  cocoa
// Anna and Michael made me check without writing a proper message.
//
// Revision 1.9  2006/11/24 17:04:32  cocoa
// -- reorganized includes of header files
//
// Revision 1.8  2006/11/23 17:39:11  cocoa
// -- added #include
//
// Revision 1.7  2006/11/16 11:27:20  cocoa
// -- reinserted myRefCountZero(): sometimes really necessary, in general safe
//
// Revision 1.6  2006/11/14 17:29:20  cocoa
// -- commented out myRefCountZero() (not necessary???)
//
// Revision 1.5  2006/10/16 23:18:59  cocoa
// Corrected use of std::swap and various special swap functions.
// Improved myApply memfn for homs of RingDistrMPolyInlPP.
//
// Revision 1.4  2006/10/06 14:04:14  cocoa
// Corrected position of #ifndef in header files.
// Separated CoCoA_ASSERT into assert.H from config.H;
// many minor consequential changes (have to #include assert.H).
// A little tidying of #include directives (esp. in Max's code).
//
// Revision 1.3  2006/08/07 21:23:25  cocoa
// Removed almost all publicly visible references to SmallExponent_t;
// changed to long in all PPMonoid functions and SparsePolyRing functions.
// DivMask remains to sorted out.
//
// Revision 1.2  2006/06/21 17:07:10  cocoa
// Fixed IsIndet bug -- why are there three almost identical copies of code?
//
// Revision 1.1.1.1  2006/05/30 11:39:37  cocoa
// Imported files
//
// Revision 1.9  2006/04/27 13:45:30  cocoa
// Changed name of NewIdentityRingHom to NewIdentityHom.
// Changed name of member functions which print out their own object
// into myOutputSelf (to distinguish from "transitive" myOutput fns).
//
// Revision 1.8  2006/03/27 12:21:25  cocoa
// Minor silly changes to reduce number of complaints from some compiler or other.
//
// Revision 1.7  2006/03/15 18:09:31  cocoa
// Changed names of member functions which print out their object
// into myOutputSelf -- hope this will appease the Intel C++ compiler.
//
// Revision 1.6  2006/03/14 17:21:18  cocoa
// Moved concrete PPMonoid impls entirely into their respective .C files.
// Now the corresponding .H files are very compact.
//
// Revision 1.5  2006/03/12 21:28:33  cocoa
// Major check in after many changes
//
// Revision 1.4  2006/02/20 22:41:20  cocoa
// All forms of the log function for power products now return SmallExponent_t
// (instead of int).  exponents now resizes the vector rather than requiring
// the user to pass in the correct size.
//
// Revision 1.3  2006/01/17 10:23:08  cocoa
// Updated DivMask; many consequential changes.
// A few other minor fixes.
//
// Revision 1.2  2005/11/17 13:25:13  cocoa
// -- "unsigned int" --> "SmallExponent_t" in ctor PPMonoidEvImpl(o, IndetNames)
//
// Revision 1.1.1.1  2005/10/17 10:46:54  cocoa
// Imported files
//
// Revision 1.9  2005/10/11 16:37:30  cocoa
// Added new small prime finite field class (see RingFpDouble).
//
// Cleaned makefiles and configuration script.
//
// Tidied PPMonoid code (to eliminate compiler warnings).
//
// Fixed bug in RingFloat::myIsInteger.
//
// Revision 1.8  2005/08/08 16:36:32  cocoa
// Just checking in before going on holiday.
// Don't really recall what changes have been made.
// Added IsIndet function for RingElem, PPMonoidElem,
// and a member function of OrdvArith.
// Improved the way failed assertions are handled.
//
// Revision 1.7  2005/07/19 15:30:20  cocoa
// A first attempt at iterators over sparse polynomials.
// Main additions are to SparsePolyRing, DistrMPoly*.
// Some consequential changes to PPMonoid*.
//
// Revision 1.6  2005/07/08 15:09:28  cocoa
// Added new symbol class (to represent names of indets).
// Integrated the new class into concrete polynomial rings
// and PPMonoid -- many consequential changes.
// Change ctors for the "inline" sparse poly rings: they no
// longer expect a PPMonoid, but build their own instead
// (has to be a PPMonoidOv).
//
// Revision 1.5  2005/07/01 16:08:15  cocoa
// Friday check-in.  Major change to structure under PolyRing:
// now SparsePolyRing and DUPolyRing are separated (in preparation
// for implementing iterators).
//
// A number of other relatively minor changes had to be chased through
// (e.g. IndetPower).
//
// Revision 1.4  2005/06/23 15:42:41  cocoa
// Fixed typo in GNU fdl -- all doc/*.txt files affected.
// Minor corrections to PPMonoid (discovered while writing doc).
//
// Revision 1.3  2005/06/22 14:47:56  cocoa
// PPMonoids and PPMonoidElems updated to mirror the structure
// used for rings and RingElems.  Many consequential changes.
//
// Revision 1.2  2005/05/04 16:50:31  cocoa
// -- new code for MatrixOrderingImpl
// -- changed: CmpBase also requires GradingDim
//
// Revision 1.1.1.1  2005/05/03 15:47:31  cocoa
// Imported files
//
// Revision 1.4  2005/04/29 15:42:02  cocoa
// Improved documentation for GMPAllocator.
// Added example program for GMPAllocator.
// Added example program for simple ops on polynomials.
// Added two new ctors for (principal) ideals (from long, and from ZZ).
// Added (crude) printing for PPMonoids.
// Updated library.H (#included GMPAllocator.H).
//
// Revision 1.3  2005/04/20 15:40:48  cocoa
// Major change: modified the standard way errors are to be signalled
// (now via a macro which records filename and line number).  Updated
// documentation in error.txt accordingly.
//
// Improved the documentation in matrix.txt (still more work to be done).
//
// Revision 1.2  2005/04/19 14:06:04  cocoa
// Added GPL and GFDL licence stuff.
//
// Revision 1.1.1.1  2005/01/27 15:12:13  cocoa
// Imported files
//
// Revision 1.3  2004/11/25 16:14:21  cocoa
// (1) Fixed definition of specialization of std::swap template function
//     so that it compiles with gcc 3.4.3
// (2) Implemented monomial function for polynomial rings.
// (3) Added one(PPM) and PPM->myOne() functions.
//
// Revision 1.2  2004/11/12 15:49:29  cocoa
// Tidying prior to 0.90 release.
// (a) documentation improved (or marked as poor)
// (b) sundry minor improvements to the code
//
// Revision 1.1  2004/11/11 13:45:22  cocoa
// -- renamed: was PPMonoidSafe
//
// Revision 1.3  2004/11/03 17:52:05  cocoa
// -- added just a warning to remind me to implement matrix ordering
//
// Revision 1.2  2004/11/02 14:56:33  cocoa
// -- changed *Print* into *Output* (myPrint --> myOutput)
// -- changed *Var* into *Indet* (myPrintVarName --> myOutputIndetName)
// -- removed suffix "IgnoreDivMask"
// -- added myComputeDivMask
// -- improved storing of IndetNames
// -- changed ExpvElem into SmallExponent_t
//
// Revision 1.1  2004/10/29 15:31:25  cocoa
// -- new PPMonoid for compatibility with OrdvArith (without DivMask)
//