X-Git-Url: https://scm.cri.ensmp.fr/git/Faustine.git/blobdiff_plain/c7f552fd8888da2f0d8cfb228fe0f28d3df3a12c..b4b6f2ea75b9f0f3ca918f5b84016610bf7a4d4f:/interpretor/preprocessor/faust-0.9.47mr3/compiler/normalize/mterm.cpp

diff --git a/interpretor/preprocessor/faust-0.9.47mr3/compiler/normalize/mterm.cpp b/interpretor/preprocessor/faust-0.9.47mr3/compiler/normalize/mterm.cpp
new file mode 100644
index 0000000..1f8bd2d
--- /dev/null
+++ b/interpretor/preprocessor/faust-0.9.47mr3/compiler/normalize/mterm.cpp
@@ -0,0 +1,485 @@
+#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) {}
+
+/**
+ * create a mterm from a tree sexpression
+ */
+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
+ */
+bool mterm::isNotZero() const
+{
+	return !isZero(fCoef);
+}
+
+/**
+ * true if mterm doesn't represent number 0
+ */
+bool mterm::isNegative() const
+{
+	return !isGEZero(fCoef);
+}
+
+/**
+ * print a mterm in a human readable format
+ */
+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)
+ */
+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
+ */
+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
+ */
+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
+ */
+const mterm& mterm::operator = (const mterm& m)
+{
+	fCoef = m.fCoef;
+	fFactors = m.fFactors;
+    return *this;
+}
+
+/**
+ * Clean a mterm by removing x**0 factors 
+ */
+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
+ */
+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
+ */
+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
+ */
+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
+ */
+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
+ */
+mterm mterm::operator * (const mterm& m) const
+{
+    mterm r(*this);
+    r *= m;
+    return r;
+}
+
+/**
+ * Divide two mterms
+ */
+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
+ */
+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
+ */
+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;
+    else exit(1);
+}
+
+/**
+ * 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);
+    else exit(1);
+}
+
+/**	
+ * 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) {
+        assert(q != 0);
+		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)
+ */
+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
+ */
+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;
+	}
+}
+