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