mterm.cpp

Go to the documentation of this file.

#include "mterm.hh"
#include "signals.hh"
#include "ppsig.hh"
#include "xtended.hh"
#include <assert.h>
//static void collectMulTerms (Tree& coef, map<Tree,int>& M, Tree t, bool invflag=false);

#undef TRACE

using namespace std;

typedef map<Tree,int> MP;

mterm::mterm ()                 : fCoef(sigInt(0)) {}
mterm::mterm (int k)            : fCoef(sigInt(k)) {}
mterm::mterm (double k)         : fCoef(sigReal(k)) {}  // cerr << "DOUBLE " << endl; }
mterm::mterm (const mterm& m)   : fCoef(m.fCoef), fFactors(m.fFactors) {}

mterm::mterm (Tree t) : fCoef(sigInt(1))
{
    //cerr << "mterm::mterm (Tree t) : " << ppsig(t) << endl;
    *this *= t; 
    //cerr << "MTERM(" << ppsig(t) << ") -> " << *this << endl;
}

bool mterm::isNotZero() const
{
    return !isZero(fCoef);
}

bool mterm::isNegative() const
{
    return !isGEZero(fCoef);
}

ostream& mterm::print(ostream& dst) const
{
    const char* sep = "";
    if (!isOne(fCoef) || fFactors.empty()) { dst << ppsig(fCoef); sep = " * "; }
    //if (true) { dst << ppsig(fCoef); sep = " * "; }
    for (MP::const_iterator p = fFactors.begin(); p != fFactors.end(); p++) {
        dst << sep << ppsig(p->first);
        if (p->second != 1) dst << "**" << p->second;
        sep = " * ";
    }
    return dst;
}

int mterm::complexity() const
{
    int c = isOne(fCoef) ? 0 : 1;
    for (MP::const_iterator p = fFactors.begin(); p != fFactors.end(); ++p) {
        c += (1+getSigOrder(p->first))*abs(p->second);
    }
    return c;
}

static bool isSigPow(Tree sig, Tree& x, int& n)
{
    //cerr << "isSigPow("<< *sig << ')' << endl;
    xtended* p = (xtended*) getUserData(sig);
    if (p == gPowPrim) {
       if (isSigInt(sig->branch(1), &n)) {
            x = sig->branch(0);
            //cerr << "factor of isSigPow " << *x << endl;
            return true;
        }
    }
    return false;
}

static Tree sigPow(Tree x, int p)
{
    return tree(gPowPrim->symbol(), x, sigInt(p));
}

const mterm& mterm::operator *= (Tree t)
{
    int     op, n;
    Tree    x,y;

    assert(t!=0);

    if (isNum(t)) {
        fCoef = mulNums(fCoef,t);

    } else if (isSigBinOp(t, &op, x, y) && (op == kMul)) {
        *this *= x;
        *this *= y;

    } else if (isSigBinOp(t, &op, x, y) && (op == kDiv)) {
        *this *= x;
        *this /= y;

    } else {
        if (isSigPow(t,x,n)) {
            fFactors[x] += n;
        } else {
            fFactors[t] += 1;
        }
    }
    return *this;
}

const mterm& mterm::operator /= (Tree t)
{
    //cerr << "division en place : " << *this << " / " << ppsig(t) << endl;
    int     op,n;
    Tree    x,y;

    assert(t!=0);

    if (isNum(t)) {
        fCoef = divExtendedNums(fCoef,t);

    } else if (isSigBinOp(t, &op, x, y) && (op == kMul)) {
        *this /= x;
        *this /= y;

    } else if (isSigBinOp(t, &op, x, y) && (op == kDiv)) {
        *this /= x;
        *this *= y;

    } else {
        if (isSigPow(t,x,n)) {
            fFactors[x] -= n;
        } else {
            fFactors[t] -= 1;
        }
    }
    return *this;
}

const mterm& mterm::operator = (const mterm& m)
{
    fCoef = m.fCoef;
    fFactors = m.fFactors;
    return *this;
}

void mterm::cleanup()
{
    if (isZero(fCoef)) {
        fFactors.clear();
    } else {
        for (MP::iterator p = fFactors.begin(); p != fFactors.end(); ) {
            if (p->second == 0) {
                fFactors.erase(p++);
            } else {
                p++;
            }
        }
    }
}

const mterm& mterm::operator += (const mterm& m)
{
    if (isZero(m.fCoef)) {
        // nothing to do
    } else if (isZero(fCoef))   {
        // copy of m
        fCoef = m.fCoef;
        fFactors = m.fFactors;
    } else {
        // only add mterms of same signature
        assert(signatureTree() == m.signatureTree());
        fCoef = addNums(fCoef, m.fCoef);
    }
    cleanup();
    return *this;
}

const mterm& mterm::operator -= (const mterm& m)
{
    if (isZero(m.fCoef)) {
        // nothing to do
    } else if (isZero(fCoef))   {
        // minus of m
        fCoef = minusNum(m.fCoef);
        fFactors = m.fFactors;
    } else {
        // only add mterms of same signature
        assert(signatureTree() == m.signatureTree());
        fCoef = subNums(fCoef, m.fCoef);
    }
    cleanup();
    return *this;
}

const mterm& mterm::operator *= (const mterm& m)
{
    fCoef = mulNums(fCoef,m.fCoef);
    for (MP::const_iterator p = m.fFactors.begin(); p != m.fFactors.end(); p++) {
        fFactors[p->first] += p->second;
    }
    cleanup();
    return *this;
}

const mterm& mterm::operator /= (const mterm& m)
{
    //cerr << "division en place : " << *this << " / " << m << endl;
    fCoef = divExtendedNums(fCoef,m.fCoef);
    for (MP::const_iterator p = m.fFactors.begin(); p != m.fFactors.end(); p++) {
        fFactors[p->first] -= p->second;
    }
    cleanup();
    return *this;
}

mterm mterm::operator * (const mterm& m) const
{
    mterm r(*this);
    r *= m;
    return r;
}

mterm mterm::operator / (const mterm& m) const
{
    mterm r(*this);
    r /= m;
    return r;
}

static int common(int a, int b)
{
    if (a > 0 & b > 0) {
        return min(a,b);
    } else if (a < 0 & b < 0) {
        return max(a,b);
    } else {
        return 0;
    }
}


mterm gcd (const mterm& m1, const mterm& m2)
{
    //cerr << "GCD of " << m1 << " and " << m2 << endl;

    Tree c = (m1.fCoef == m2.fCoef) ? m1.fCoef : tree(1);       // common coefficient (real gcd not needed)
    mterm R(c);
    for (MP::const_iterator p1 = m1.fFactors.begin(); p1 != m1.fFactors.end(); p1++) {
        Tree t = p1->first;
        MP::const_iterator p2 = m2.fFactors.find(t);
        if (p2 != m2.fFactors.end()) {
            int v1 = p1->second;
            int v2 = p2->second;
            int c = common(v1,v2);
            if (c != 0) {
                R.fFactors[t] = c;
            }
        }
    }
    //cerr << "GCD of " << m1 << " and " << m2 << " is : " << R << endl;
    return R;
}

static bool contains(int a, int b)
{
    return (b == 0) || (a/b > 0);
}

bool mterm::hasDivisor (const mterm& n) const
{
    for (MP::const_iterator p1 = n.fFactors.begin(); p1 != n.fFactors.end(); p1++) {
        // for each factor f**q of m
        Tree    f = p1->first;  
        int     v = p1->second;
        
        // check that f is also a factor of *this
        MP::const_iterator p2 = fFactors.find(f);
        if (p2 == fFactors.end()) return false;
        
        // analyze quantities
        int u = p2->second;
        if (! contains(u,v) ) return false;
    }
    return true;
}

static Tree buildPowTerm(Tree f, int q)
{
    assert(f);
    assert(q>0);
    if (q>1) {
        return sigPow(f, q);
    } else {
        return f;
    }
}

static void combineMulLeft(Tree& R, Tree A)
{
    if (R && A)     R = sigMul(R,A);
    else if (A)     R = A;
    else exit(1);
}

static void combineDivLeft(Tree& R, Tree A)
{
    if (R && A)     R = sigDiv(R,A);
    else if (A)     R = sigDiv(tree(1.0f),A);
    else exit(1);
}

static void combineMulDiv(Tree& M, Tree& D, Tree f, int q)
{
    #ifdef TRACE
    cerr << "combineMulDiv (" << M << "/"  << D << "*" << ppsig(f)<< "**" << q << endl;
    #endif
    if (f) {
        assert(q != 0);
        if (q > 0) {
            combineMulLeft(M, buildPowTerm(f,q));
        } else if (q < 0) {
            combineMulLeft(D, buildPowTerm(f,-q));
        }
    }
}   
    
            
Tree mterm::signatureTree() const
{
    return normalizedTree(true);
}
    
Tree mterm::normalizedTree(bool signatureMode, bool negativeMode) const
{
    if (fFactors.empty() || isZero(fCoef)) {
        // it's a pure number
        if (signatureMode)  return tree(1);
        if (negativeMode)   return minusNum(fCoef);
        else                return fCoef;
    } else {
        // it's not a pure number, it has factors
        Tree A[4], B[4];
        
        // group by order
        for (int order = 0; order < 4; order++) {
            A[order] = 0; B[order] = 0;
            for (MP::const_iterator p = fFactors.begin(); p != fFactors.end(); p++) {
                Tree    f = p->first;       // f = factor
                int     q = p->second;      // q = power of f
                if (f && q && getSigOrder(f)==order) {
                    
                    combineMulDiv (A[order], B[order], f, q);
                }
            }
        }
        if (A[0] != 0) cerr << "A[0] == " << *A[0] << endl; 
        if (B[0] != 0) cerr << "B[0] == " << *B[0] << endl; 
        // en principe ici l'order zero est vide car il correspond au coef numerique
        assert(A[0] == 0);
        assert(B[0] == 0);
        
        // we only use a coeficient if it differes from 1 and if we are not in signature mode
        if (! (signatureMode | isOne(fCoef))) {
            A[0] = (negativeMode) ? minusNum(fCoef) : fCoef;
        }
        
        if (signatureMode) {
            A[0] = 0;
        } else if (negativeMode) {
            if (isMinusOne(fCoef)) { A[0] = 0; } else { A[0] = minusNum(fCoef); }
        } else if (isOne(fCoef)) {
            A[0] = 0;
        } else {
            A[0] = fCoef;
        }
                    
        // combine each order separately : R[i] = A[i]/B[i]
        Tree RR = 0;
        for (int order = 0; order < 4; order++) {
            if (A[order] && B[order])   combineMulLeft(RR,sigDiv(A[order],B[order]));
            else if (A[order])          combineMulLeft(RR,A[order]);
            else if (B[order])          combineDivLeft(RR,B[order]);
        }
        if (RR == 0) RR = tree(1); // a verifier *******************
            
        assert(RR);
        //cerr << "Normalized Tree of " << *this << " is " << ppsig(RR) << endl;
        return RR;
    }
}