00001 #include "mterm.hh"
00002 #include "signals.hh"
00003 #include "ppsig.hh"
00004 #include "xtended.hh"
00005 #include <assert.h>
00006
00007
00008 #undef TRACE
00009
00010 using namespace std;
00011
00012 typedef map<Tree,int> MP;
00013
00014 mterm::mterm () : fCoef(sigInt(0)) {}
00015 mterm::mterm (int k) : fCoef(sigInt(k)) {}
00016 mterm::mterm (double k) : fCoef(sigReal(k)) {}
00017 mterm::mterm (const mterm& m) : fCoef(m.fCoef), fFactors(m.fFactors) {}
00018
00022 mterm::mterm (Tree t) : fCoef(sigInt(1))
00023 {
00024
00025 *this *= t;
00026
00027 }
00028
00032 bool mterm::isNotZero() const
00033 {
00034 return !isZero(fCoef);
00035 }
00036
00040 bool mterm::isNegative() const
00041 {
00042 return !isGEZero(fCoef);
00043 }
00044
00048 ostream& mterm::print(ostream& dst) const
00049 {
00050 const char* sep = "";
00051 if (!isOne(fCoef) || fFactors.empty()) { dst << ppsig(fCoef); sep = " * "; }
00052
00053 for (MP::const_iterator p = fFactors.begin(); p != fFactors.end(); p++) {
00054 dst << sep << ppsig(p->first);
00055 if (p->second != 1) dst << "**" << p->second;
00056 sep = " * ";
00057 }
00058 return dst;
00059 }
00060
00065 int mterm::complexity() const
00066 {
00067 int c = isOne(fCoef) ? 0 : 1;
00068 for (MP::const_iterator p = fFactors.begin(); p != fFactors.end(); ++p) {
00069 c += (1+getSigOrder(p->first))*abs(p->second);
00070 }
00071 return c;
00072 }
00073
00077 static bool isSigPow(Tree sig, Tree& x, int& n)
00078 {
00079
00080 xtended* p = (xtended*) getUserData(sig);
00081 if (p == gPowPrim) {
00082 if (isSigInt(sig->branch(1), &n)) {
00083 x = sig->branch(0);
00084
00085 return true;
00086 }
00087 }
00088 return false;
00089 }
00090
00094 static Tree sigPow(Tree x, int p)
00095 {
00096 return tree(gPowPrim->symbol(), x, sigInt(p));
00097 }
00098
00103 const mterm& mterm::operator *= (Tree t)
00104 {
00105 int op, n;
00106 Tree x,y;
00107
00108 assert(t!=0);
00109
00110 if (isNum(t)) {
00111 fCoef = mulNums(fCoef,t);
00112
00113 } else if (isSigBinOp(t, &op, x, y) && (op == kMul)) {
00114 *this *= x;
00115 *this *= y;
00116
00117 } else if (isSigBinOp(t, &op, x, y) && (op == kDiv)) {
00118 *this *= x;
00119 *this /= y;
00120
00121 } else {
00122 if (isSigPow(t,x,n)) {
00123 fFactors[x] += n;
00124 } else {
00125 fFactors[t] += 1;
00126 }
00127 }
00128 return *this;
00129 }
00130
00135 const mterm& mterm::operator /= (Tree t)
00136 {
00137
00138 int op,n;
00139 Tree x,y;
00140
00141 assert(t!=0);
00142
00143 if (isNum(t)) {
00144 fCoef = divExtendedNums(fCoef,t);
00145
00146 } else if (isSigBinOp(t, &op, x, y) && (op == kMul)) {
00147 *this /= x;
00148 *this /= y;
00149
00150 } else if (isSigBinOp(t, &op, x, y) && (op == kDiv)) {
00151 *this /= x;
00152 *this *= y;
00153
00154 } else {
00155 if (isSigPow(t,x,n)) {
00156 fFactors[x] -= n;
00157 } else {
00158 fFactors[t] -= 1;
00159 }
00160 }
00161 return *this;
00162 }
00163
00167 const mterm& mterm::operator = (const mterm& m)
00168 {
00169 fCoef = m.fCoef;
00170 fFactors = m.fFactors;
00171 return *this;
00172 }
00173
00177 void mterm::cleanup()
00178 {
00179 if (isZero(fCoef)) {
00180 fFactors.clear();
00181 } else {
00182 for (MP::iterator p = fFactors.begin(); p != fFactors.end(); ) {
00183 if (p->second == 0) {
00184 fFactors.erase(p++);
00185 } else {
00186 p++;
00187 }
00188 }
00189 }
00190 }
00191
00196 const mterm& mterm::operator += (const mterm& m)
00197 {
00198 if (isZero(m.fCoef)) {
00199
00200 } else if (isZero(fCoef)) {
00201
00202 fCoef = m.fCoef;
00203 fFactors = m.fFactors;
00204 } else {
00205
00206 assert(signatureTree() == m.signatureTree());
00207 fCoef = addNums(fCoef, m.fCoef);
00208 }
00209 cleanup();
00210 return *this;
00211 }
00212
00217 const mterm& mterm::operator -= (const mterm& m)
00218 {
00219 if (isZero(m.fCoef)) {
00220
00221 } else if (isZero(fCoef)) {
00222
00223 fCoef = minusNum(m.fCoef);
00224 fFactors = m.fFactors;
00225 } else {
00226
00227 assert(signatureTree() == m.signatureTree());
00228 fCoef = subNums(fCoef, m.fCoef);
00229 }
00230 cleanup();
00231 return *this;
00232 }
00233
00237 const mterm& mterm::operator *= (const mterm& m)
00238 {
00239 fCoef = mulNums(fCoef,m.fCoef);
00240 for (MP::const_iterator p = m.fFactors.begin(); p != m.fFactors.end(); p++) {
00241 fFactors[p->first] += p->second;
00242 }
00243 cleanup();
00244 return *this;
00245 }
00246
00250 const mterm& mterm::operator /= (const mterm& m)
00251 {
00252
00253 fCoef = divExtendedNums(fCoef,m.fCoef);
00254 for (MP::const_iterator p = m.fFactors.begin(); p != m.fFactors.end(); p++) {
00255 fFactors[p->first] -= p->second;
00256 }
00257 cleanup();
00258 return *this;
00259 }
00260
00264 mterm mterm::operator * (const mterm& m) const
00265 {
00266 mterm r(*this);
00267 r *= m;
00268 return r;
00269 }
00270
00274 mterm mterm::operator / (const mterm& m) const
00275 {
00276 mterm r(*this);
00277 r /= m;
00278 return r;
00279 }
00280
00284 static int common(int a, int b)
00285 {
00286 if (a > 0 & b > 0) {
00287 return min(a,b);
00288 } else if (a < 0 & b < 0) {
00289 return max(a,b);
00290 } else {
00291 return 0;
00292 }
00293 }
00294
00295
00299 mterm gcd (const mterm& m1, const mterm& m2)
00300 {
00301
00302
00303 Tree c = (m1.fCoef == m2.fCoef) ? m1.fCoef : tree(1);
00304 mterm R(c);
00305 for (MP::const_iterator p1 = m1.fFactors.begin(); p1 != m1.fFactors.end(); p1++) {
00306 Tree t = p1->first;
00307 MP::const_iterator p2 = m2.fFactors.find(t);
00308 if (p2 != m2.fFactors.end()) {
00309 int v1 = p1->second;
00310 int v2 = p2->second;
00311 int c = common(v1,v2);
00312 if (c != 0) {
00313 R.fFactors[t] = c;
00314 }
00315 }
00316 }
00317
00318 return R;
00319 }
00320
00325 static bool contains(int a, int b)
00326 {
00327 return (b == 0) || (a/b > 0);
00328 }
00329
00337 bool mterm::hasDivisor (const mterm& n) const
00338 {
00339 for (MP::const_iterator p1 = n.fFactors.begin(); p1 != n.fFactors.end(); p1++) {
00340
00341 Tree f = p1->first;
00342 int v = p1->second;
00343
00344
00345 MP::const_iterator p2 = fFactors.find(f);
00346 if (p2 == fFactors.end()) return false;
00347
00348
00349 int u = p2->second;
00350 if (! contains(u,v) ) return false;
00351 }
00352 return true;
00353 }
00354
00362 static Tree buildPowTerm(Tree f, int q)
00363 {
00364 assert(f);
00365 assert(q>0);
00366 if (q>1) {
00367 return sigPow(f, q);
00368 } else {
00369 return f;
00370 }
00371 }
00372
00376 static void combineMulLeft(Tree& R, Tree A)
00377 {
00378 if (R && A) R = sigMul(R,A);
00379 else if (A) R = A;
00380 }
00381
00385 static void combineDivLeft(Tree& R, Tree A)
00386 {
00387 if (R && A) R = sigDiv(R,A);
00388 else if (A) R = sigDiv(tree(1.0f),A);
00389 }
00390
00394 static void combineMulDiv(Tree& M, Tree& D, Tree f, int q)
00395 {
00396 #ifdef TRACE
00397 cerr << "combineMulDiv (" << M << "/" << D << "*" << ppsig(f)<< "**" << q << endl;
00398 #endif
00399 if (f) {
00400 if (q > 0) {
00401 combineMulLeft(M, buildPowTerm(f,q));
00402 } else if (q < 0) {
00403 combineMulLeft(D, buildPowTerm(f,-q));
00404 }
00405 }
00406 }
00407
00408
00413 Tree mterm::signatureTree() const
00414 {
00415 return normalizedTree(true);
00416 }
00417
00424 Tree mterm::normalizedTree(bool signatureMode, bool negativeMode) const
00425 {
00426 if (fFactors.empty() || isZero(fCoef)) {
00427
00428 if (signatureMode) return tree(1);
00429 if (negativeMode) return minusNum(fCoef);
00430 else return fCoef;
00431 } else {
00432
00433 Tree A[4], B[4];
00434
00435
00436 for (int order = 0; order < 4; order++) {
00437 A[order] = 0; B[order] = 0;
00438 for (MP::const_iterator p = fFactors.begin(); p != fFactors.end(); p++) {
00439 Tree f = p->first;
00440 int q = p->second;
00441 if (f && q && getSigOrder(f)==order) {
00442
00443 combineMulDiv (A[order], B[order], f, q);
00444 }
00445 }
00446 }
00447 if (A[0] != 0) cerr << "A[0] == " << *A[0] << endl;
00448 if (B[0] != 0) cerr << "B[0] == " << *B[0] << endl;
00449
00450 assert(A[0] == 0);
00451 assert(B[0] == 0);
00452
00453
00454 if (! (signatureMode | isOne(fCoef))) {
00455 A[0] = (negativeMode) ? minusNum(fCoef) : fCoef;
00456 }
00457
00458 if (signatureMode) {
00459 A[0] = 0;
00460 } else if (negativeMode) {
00461 if (isMinusOne(fCoef)) { A[0] = 0; } else { A[0] = minusNum(fCoef); }
00462 } else if (isOne(fCoef)) {
00463 A[0] = 0;
00464 } else {
00465 A[0] = fCoef;
00466 }
00467
00468
00469 Tree RR = 0;
00470 for (int order = 0; order < 4; order++) {
00471 if (A[order] && B[order]) combineMulLeft(RR,sigDiv(A[order],B[order]));
00472 else if (A[order]) combineMulLeft(RR,A[order]);
00473 else if (B[order]) combineDivLeft(RR,B[order]);
00474 }
00475 if (RR == 0) RR = tree(1);
00476
00477 assert(RR);
00478
00479 return RR;
00480 }
00481 }
00482