Logo Search packages:      
Sourcecode: faust version File versions  Download package

mterm.cpp

#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;

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

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

/**
 * true if mterm doesn't represent number 0
 */
00032 bool mterm::isNotZero() const
{
      return !isZero(fCoef);
}

/**
 * true if mterm doesn't represent number 0
 */
00040 bool mterm::isNegative() const
{
      return !isGEZero(fCoef);
}

/**
 * print a mterm in a human readable format
 */
00048 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;
}

/**
 * Compute the "complexity" of a mterm, that is the number of
 * factors it contains (weighted by the importance of these factors)
 */
00065 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;
}

/**
 * match x^p with p:int
 */
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;
}

/**
 * produce x^p with p:int
 */
static Tree sigPow(Tree x, int p)
{
      return tree(gPowPrim->symbol(), x, sigInt(p));
}

/**
 * Multiple a mterm by an expression tree t. Go down recursively looking 
 * for multiplications and divisions
 */
00103 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;
}

/**
 * Divide a mterm by an expression tree t. Go down recursively looking 
 * for multiplications and divisions
 */
00135 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;
}

/**
 * replace the content with a copy of m
 */
00167 const mterm& mterm::operator = (const mterm& m)
{
      fCoef = m.fCoef;
      fFactors = m.fFactors;
    return *this;
}

/**
 * Clean a mterm by removing x**0 factors 
 */
00177 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++;
                  }
            }
      }
}

/**
 * Add in place an mterm. As we want the result to be
 * a mterm therefore essentially mterms of same signature can be added
 */
00196 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;
}

/**
 * Substract in place an mterm. As we want the result to be
 * a mterm therefore essentially mterms of same signature can be substracted
 */
00217 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;
}

/**
 * Multiply a mterm by the content of another mterm
 */
00237 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;
}

/**
 * Divide a mterm by the content of another mterm
 */
00250 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;
}

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

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

/**
 * return the "common quantity" of two numbers
 */
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;
    }
}


/**
 * return a mterm that is the greatest common divisor of two mterms
 */
00299 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;
}

/**
 * We say that a "contains" b if a/b > 0. For example 3 contains 2 and
 * -4 contains -2, but 3 doesn't contains -2 and -3 doesn't contains 1
 */
static bool contains(int a, int b)
{
      return (b == 0) || (a/b > 0);
}

/**
 * Check if M accept N has a divisor. We can say that N is
 * a divisor of M if M = N*Q and the complexity is preserved :
 * complexity(M) = complexity(N)+complexity(Q)
 * x**u has divisor x**v if u >= v
 * x**-u has divisor x**-v if -u <= -v
 */
00337 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;
}

/**
 * produce the canonical tree correspoding to a mterm
 */
 
/**
 * Build a power term of type f**q -> (((f.f).f)..f) with q>0
 */
static Tree buildPowTerm(Tree f, int q)
{
      assert(f);
      assert(q>0);
      if (q>1) {
            return sigPow(f, q);
      } else {
            return f;
      }
}

/**
 * Combine R and A doing R = R*A or R = A
 */
static void combineMulLeft(Tree& R, Tree A)
{
      if (R && A)       R = sigMul(R,A);
      else if (A)       R = A;
}

/**
 * Combine R and A doing R = R*A or R = A
 */
static void combineDivLeft(Tree& R, Tree A)
{
      if (R && A)       R = sigDiv(R,A);
      else if (A)       R = sigDiv(tree(1.0f),A);
}

/**   
 * Do M = M * f**q or D = D * f**-q
 */
static void combineMulDiv(Tree& M, Tree& D, Tree f, int q)
{
      #ifdef TRACE
      cerr << "combineMulDiv (" << M << "/"  << D << "*" << ppsig(f)<< "**" << q << endl;
      #endif
      if (f) {
            if (q > 0) {
                  combineMulLeft(M, buildPowTerm(f,q));
            } else if (q < 0) {
                  combineMulLeft(D, buildPowTerm(f,-q));
            }
      }
}     
      
                  
/**
 * returns a normalized (canonical) tree expression of structure :
 *          ((v1/v2)*(c1/c2))*(s1/s2)
 */
00413 Tree mterm::signatureTree() const
{
      return normalizedTree(true);
}
      
/**
 * returns a normalized (canonical) tree expression of structure :
 *          ((k*(v1/v2))*(c1/c2))*(s1/s2)
 * In signature mode the fCoef factor is ommited
 * In negativeMode the sign of the fCoef factor is inverted
 */
00424 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;
      }
}


Generated by  Doxygen 1.6.0   Back to index