Merge branch 'master' of https://scm.cri.ensmp.fr/git/Faustine
[Faustine.git] / interpretor / faust-0.9.47mr3 / compiler / normalize / mterm.cpp
1 #include "mterm.hh"
2 #include "signals.hh"
3 #include "ppsig.hh"
4 #include "xtended.hh"
5 #include <assert.h>
6 //static void collectMulTerms (Tree& coef, map<Tree,int>& M, Tree t, bool invflag=false);
7
8 #undef TRACE
9
10 using namespace std;
11
12 typedef map<Tree,int> MP;
13
14 mterm::mterm () : fCoef(sigInt(0)) {}
15 mterm::mterm (int k) : fCoef(sigInt(k)) {}
16 mterm::mterm (double k) : fCoef(sigReal(k)) {} // cerr << "DOUBLE " << endl; }
17 mterm::mterm (const mterm& m) : fCoef(m.fCoef), fFactors(m.fFactors) {}
18
19 /**
20 * create a mterm from a tree sexpression
21 */
22 mterm::mterm (Tree t) : fCoef(sigInt(1))
23 {
24 //cerr << "mterm::mterm (Tree t) : " << ppsig(t) << endl;
25 *this *= t;
26 //cerr << "MTERM(" << ppsig(t) << ") -> " << *this << endl;
27 }
28
29 /**
30 * true if mterm doesn't represent number 0
31 */
32 bool mterm::isNotZero() const
33 {
34 return !isZero(fCoef);
35 }
36
37 /**
38 * true if mterm doesn't represent number 0
39 */
40 bool mterm::isNegative() const
41 {
42 return !isGEZero(fCoef);
43 }
44
45 /**
46 * print a mterm in a human readable format
47 */
48 ostream& mterm::print(ostream& dst) const
49 {
50 const char* sep = "";
51 if (!isOne(fCoef) || fFactors.empty()) { dst << ppsig(fCoef); sep = " * "; }
52 //if (true) { dst << ppsig(fCoef); sep = " * "; }
53 for (MP::const_iterator p = fFactors.begin(); p != fFactors.end(); p++) {
54 dst << sep << ppsig(p->first);
55 if (p->second != 1) dst << "**" << p->second;
56 sep = " * ";
57 }
58 return dst;
59 }
60
61 /**
62 * Compute the "complexity" of a mterm, that is the number of
63 * factors it contains (weighted by the importance of these factors)
64 */
65 int mterm::complexity() const
66 {
67 int c = isOne(fCoef) ? 0 : 1;
68 for (MP::const_iterator p = fFactors.begin(); p != fFactors.end(); ++p) {
69 c += (1+getSigOrder(p->first))*abs(p->second);
70 }
71 return c;
72 }
73
74 /**
75 * match x^p with p:int
76 */
77 static bool isSigPow(Tree sig, Tree& x, int& n)
78 {
79 //cerr << "isSigPow("<< *sig << ')' << endl;
80 xtended* p = (xtended*) getUserData(sig);
81 if (p == gPowPrim) {
82 if (isSigInt(sig->branch(1), &n)) {
83 x = sig->branch(0);
84 //cerr << "factor of isSigPow " << *x << endl;
85 return true;
86 }
87 }
88 return false;
89 }
90
91 /**
92 * produce x^p with p:int
93 */
94 static Tree sigPow(Tree x, int p)
95 {
96 return tree(gPowPrim->symbol(), x, sigInt(p));
97 }
98
99 /**
100 * Multiple a mterm by an expression tree t. Go down recursively looking
101 * for multiplications and divisions
102 */
103 const mterm& mterm::operator *= (Tree t)
104 {
105 int op, n;
106 Tree x,y;
107
108 assert(t!=0);
109
110 if (isNum(t)) {
111 fCoef = mulNums(fCoef,t);
112
113 } else if (isSigBinOp(t, &op, x, y) && (op == kMul)) {
114 *this *= x;
115 *this *= y;
116
117 } else if (isSigBinOp(t, &op, x, y) && (op == kDiv)) {
118 *this *= x;
119 *this /= y;
120
121 } else {
122 if (isSigPow(t,x,n)) {
123 fFactors[x] += n;
124 } else {
125 fFactors[t] += 1;
126 }
127 }
128 return *this;
129 }
130
131 /**
132 * Divide a mterm by an expression tree t. Go down recursively looking
133 * for multiplications and divisions
134 */
135 const mterm& mterm::operator /= (Tree t)
136 {
137 //cerr << "division en place : " << *this << " / " << ppsig(t) << endl;
138 int op,n;
139 Tree x,y;
140
141 assert(t!=0);
142
143 if (isNum(t)) {
144 fCoef = divExtendedNums(fCoef,t);
145
146 } else if (isSigBinOp(t, &op, x, y) && (op == kMul)) {
147 *this /= x;
148 *this /= y;
149
150 } else if (isSigBinOp(t, &op, x, y) && (op == kDiv)) {
151 *this /= x;
152 *this *= y;
153
154 } else {
155 if (isSigPow(t,x,n)) {
156 fFactors[x] -= n;
157 } else {
158 fFactors[t] -= 1;
159 }
160 }
161 return *this;
162 }
163
164 /**
165 * replace the content with a copy of m
166 */
167 const mterm& mterm::operator = (const mterm& m)
168 {
169 fCoef = m.fCoef;
170 fFactors = m.fFactors;
171 return *this;
172 }
173
174 /**
175 * Clean a mterm by removing x**0 factors
176 */
177 void mterm::cleanup()
178 {
179 if (isZero(fCoef)) {
180 fFactors.clear();
181 } else {
182 for (MP::iterator p = fFactors.begin(); p != fFactors.end(); ) {
183 if (p->second == 0) {
184 fFactors.erase(p++);
185 } else {
186 p++;
187 }
188 }
189 }
190 }
191
192 /**
193 * Add in place an mterm. As we want the result to be
194 * a mterm therefore essentially mterms of same signature can be added
195 */
196 const mterm& mterm::operator += (const mterm& m)
197 {
198 if (isZero(m.fCoef)) {
199 // nothing to do
200 } else if (isZero(fCoef)) {
201 // copy of m
202 fCoef = m.fCoef;
203 fFactors = m.fFactors;
204 } else {
205 // only add mterms of same signature
206 assert(signatureTree() == m.signatureTree());
207 fCoef = addNums(fCoef, m.fCoef);
208 }
209 cleanup();
210 return *this;
211 }
212
213 /**
214 * Substract in place an mterm. As we want the result to be
215 * a mterm therefore essentially mterms of same signature can be substracted
216 */
217 const mterm& mterm::operator -= (const mterm& m)
218 {
219 if (isZero(m.fCoef)) {
220 // nothing to do
221 } else if (isZero(fCoef)) {
222 // minus of m
223 fCoef = minusNum(m.fCoef);
224 fFactors = m.fFactors;
225 } else {
226 // only add mterms of same signature
227 assert(signatureTree() == m.signatureTree());
228 fCoef = subNums(fCoef, m.fCoef);
229 }
230 cleanup();
231 return *this;
232 }
233
234 /**
235 * Multiply a mterm by the content of another mterm
236 */
237 const mterm& mterm::operator *= (const mterm& m)
238 {
239 fCoef = mulNums(fCoef,m.fCoef);
240 for (MP::const_iterator p = m.fFactors.begin(); p != m.fFactors.end(); p++) {
241 fFactors[p->first] += p->second;
242 }
243 cleanup();
244 return *this;
245 }
246
247 /**
248 * Divide a mterm by the content of another mterm
249 */
250 const mterm& mterm::operator /= (const mterm& m)
251 {
252 //cerr << "division en place : " << *this << " / " << m << endl;
253 fCoef = divExtendedNums(fCoef,m.fCoef);
254 for (MP::const_iterator p = m.fFactors.begin(); p != m.fFactors.end(); p++) {
255 fFactors[p->first] -= p->second;
256 }
257 cleanup();
258 return *this;
259 }
260
261 /**
262 * Multiply two mterms
263 */
264 mterm mterm::operator * (const mterm& m) const
265 {
266 mterm r(*this);
267 r *= m;
268 return r;
269 }
270
271 /**
272 * Divide two mterms
273 */
274 mterm mterm::operator / (const mterm& m) const
275 {
276 mterm r(*this);
277 r /= m;
278 return r;
279 }
280
281 /**
282 * return the "common quantity" of two numbers
283 */
284 static int common(int a, int b)
285 {
286 if (a > 0 & b > 0) {
287 return min(a,b);
288 } else if (a < 0 & b < 0) {
289 return max(a,b);
290 } else {
291 return 0;
292 }
293 }
294
295
296 /**
297 * return a mterm that is the greatest common divisor of two mterms
298 */
299 mterm gcd (const mterm& m1, const mterm& m2)
300 {
301 //cerr << "GCD of " << m1 << " and " << m2 << endl;
302
303 Tree c = (m1.fCoef == m2.fCoef) ? m1.fCoef : tree(1); // common coefficient (real gcd not needed)
304 mterm R(c);
305 for (MP::const_iterator p1 = m1.fFactors.begin(); p1 != m1.fFactors.end(); p1++) {
306 Tree t = p1->first;
307 MP::const_iterator p2 = m2.fFactors.find(t);
308 if (p2 != m2.fFactors.end()) {
309 int v1 = p1->second;
310 int v2 = p2->second;
311 int c = common(v1,v2);
312 if (c != 0) {
313 R.fFactors[t] = c;
314 }
315 }
316 }
317 //cerr << "GCD of " << m1 << " and " << m2 << " is : " << R << endl;
318 return R;
319 }
320
321 /**
322 * We say that a "contains" b if a/b > 0. For example 3 contains 2 and
323 * -4 contains -2, but 3 doesn't contains -2 and -3 doesn't contains 1
324 */
325 static bool contains(int a, int b)
326 {
327 return (b == 0) || (a/b > 0);
328 }
329
330 /**
331 * Check if M accept N has a divisor. We can say that N is
332 * a divisor of M if M = N*Q and the complexity is preserved :
333 * complexity(M) = complexity(N)+complexity(Q)
334 * x**u has divisor x**v if u >= v
335 * x**-u has divisor x**-v if -u <= -v
336 */
337 bool mterm::hasDivisor (const mterm& n) const
338 {
339 for (MP::const_iterator p1 = n.fFactors.begin(); p1 != n.fFactors.end(); p1++) {
340 // for each factor f**q of m
341 Tree f = p1->first;
342 int v = p1->second;
343
344 // check that f is also a factor of *this
345 MP::const_iterator p2 = fFactors.find(f);
346 if (p2 == fFactors.end()) return false;
347
348 // analyze quantities
349 int u = p2->second;
350 if (! contains(u,v) ) return false;
351 }
352 return true;
353 }
354
355 /**
356 * produce the canonical tree correspoding to a mterm
357 */
358
359 /**
360 * Build a power term of type f**q -> (((f.f).f)..f) with q>0
361 */
362 static Tree buildPowTerm(Tree f, int q)
363 {
364 assert(f);
365 assert(q>0);
366 if (q>1) {
367 return sigPow(f, q);
368 } else {
369 return f;
370 }
371 }
372
373 /**
374 * Combine R and A doing R = R*A or R = A
375 */
376 static void combineMulLeft(Tree& R, Tree A)
377 {
378 if (R && A) R = sigMul(R,A);
379 else if (A) R = A;
380 else exit(1);
381 }
382
383 /**
384 * Combine R and A doing R = R*A or R = A
385 */
386 static void combineDivLeft(Tree& R, Tree A)
387 {
388 if (R && A) R = sigDiv(R,A);
389 else if (A) R = sigDiv(tree(1.0f),A);
390 else exit(1);
391 }
392
393 /**
394 * Do M = M * f**q or D = D * f**-q
395 */
396 static void combineMulDiv(Tree& M, Tree& D, Tree f, int q)
397 {
398 #ifdef TRACE
399 cerr << "combineMulDiv (" << M << "/" << D << "*" << ppsig(f)<< "**" << q << endl;
400 #endif
401 if (f) {
402 assert(q != 0);
403 if (q > 0) {
404 combineMulLeft(M, buildPowTerm(f,q));
405 } else if (q < 0) {
406 combineMulLeft(D, buildPowTerm(f,-q));
407 }
408 }
409 }
410
411
412 /**
413 * returns a normalized (canonical) tree expression of structure :
414 * ((v1/v2)*(c1/c2))*(s1/s2)
415 */
416 Tree mterm::signatureTree() const
417 {
418 return normalizedTree(true);
419 }
420
421 /**
422 * returns a normalized (canonical) tree expression of structure :
423 * ((k*(v1/v2))*(c1/c2))*(s1/s2)
424 * In signature mode the fCoef factor is ommited
425 * In negativeMode the sign of the fCoef factor is inverted
426 */
427 Tree mterm::normalizedTree(bool signatureMode, bool negativeMode) const
428 {
429 if (fFactors.empty() || isZero(fCoef)) {
430 // it's a pure number
431 if (signatureMode) return tree(1);
432 if (negativeMode) return minusNum(fCoef);
433 else return fCoef;
434 } else {
435 // it's not a pure number, it has factors
436 Tree A[4], B[4];
437
438 // group by order
439 for (int order = 0; order < 4; order++) {
440 A[order] = 0; B[order] = 0;
441 for (MP::const_iterator p = fFactors.begin(); p != fFactors.end(); p++) {
442 Tree f = p->first; // f = factor
443 int q = p->second; // q = power of f
444 if (f && q && getSigOrder(f)==order) {
445
446 combineMulDiv (A[order], B[order], f, q);
447 }
448 }
449 }
450 if (A[0] != 0) cerr << "A[0] == " << *A[0] << endl;
451 if (B[0] != 0) cerr << "B[0] == " << *B[0] << endl;
452 // en principe ici l'order zero est vide car il correspond au coef numerique
453 assert(A[0] == 0);
454 assert(B[0] == 0);
455
456 // we only use a coeficient if it differes from 1 and if we are not in signature mode
457 if (! (signatureMode | isOne(fCoef))) {
458 A[0] = (negativeMode) ? minusNum(fCoef) : fCoef;
459 }
460
461 if (signatureMode) {
462 A[0] = 0;
463 } else if (negativeMode) {
464 if (isMinusOne(fCoef)) { A[0] = 0; } else { A[0] = minusNum(fCoef); }
465 } else if (isOne(fCoef)) {
466 A[0] = 0;
467 } else {
468 A[0] = fCoef;
469 }
470
471 // combine each order separately : R[i] = A[i]/B[i]
472 Tree RR = 0;
473 for (int order = 0; order < 4; order++) {
474 if (A[order] && B[order]) combineMulLeft(RR,sigDiv(A[order],B[order]));
475 else if (A[order]) combineMulLeft(RR,A[order]);
476 else if (B[order]) combineDivLeft(RR,B[order]);
477 }
478 if (RR == 0) RR = tree(1); // a verifier *******************
479
480 assert(RR);
481 //cerr << "Normalized Tree of " << *this << " is " << ppsig(RR) << endl;
482 return RR;
483 }
484 }
485