From: Richard Kreckel Date: Sat, 29 Apr 2000 02:15:55 +0000 (+0000) Subject: - introduced info_flag::algebraic. X-Git-Tag: release_0-6-0~10 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=e009cb7984c4971df3b9e036eaead640095f46d5 - introduced info_flag::algebraic. - matrix::fraction_free_elimination(bool) finished and used by matrix::determinant(). - introduced some comments. --- diff --git a/ginac/add.cpp b/ginac/add.cpp index 717279b5..8a7e9041 100644 --- a/ginac/add.cpp +++ b/ginac/add.cpp @@ -290,21 +290,28 @@ void add::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) cons bool add::info(unsigned inf) const { - // TODO: optimize - if (inf==info_flags::polynomial || - inf==info_flags::integer_polynomial || - inf==info_flags::cinteger_polynomial || - inf==info_flags::rational_polynomial || - inf==info_flags::crational_polynomial || - inf==info_flags::rational_function) { - for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) { - if (!(recombine_pair_to_ex(*it).info(inf))) - return false; + switch (inf) { + case info_flags::polynomial: + case info_flags::integer_polynomial: + case info_flags::cinteger_polynomial: + case info_flags::rational_polynomial: + case info_flags::crational_polynomial: + case info_flags::rational_function: { + for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) { + if (!(recombine_pair_to_ex(*i).info(inf))) + return false; + } + return overall_coeff.info(inf); + } + case info_flags::algebraic: { + for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) { + if ((recombine_pair_to_ex(*i).info(inf))) + return true; + } + return false; } - return overall_coeff.info(inf); - } else { - return inherited::info(inf); } + return inherited::info(inf); } int add::degree(const symbol & s) const diff --git a/ginac/basic.cpp b/ginac/basic.cpp index 2a4a3e39..83b777bb 100644 --- a/ginac/basic.cpp +++ b/ginac/basic.cpp @@ -284,12 +284,14 @@ ex basic::collect(const symbol & s) const /* Perform automatic symbolic evaluations on expression. */ ex basic::eval(int level) const { + // There is nothing to do for basic objects: return this->hold(); } /** Evaluate object numerically. */ ex basic::evalf(int level) const { + // There is nothing to do for basic objects: return *this; } diff --git a/ginac/flags.h b/ginac/flags.h index 16228a2d..8b9caed4 100644 --- a/ginac/flags.h +++ b/ginac/flags.h @@ -45,7 +45,7 @@ public: /** Possible attributes an object can have. */ class info_flags { public: - enum { + enum { // answered by class numeric numeric, real, @@ -62,7 +62,7 @@ public: even, odd, prime, - + // answered by class relation relation, relation_equal, @@ -71,16 +71,16 @@ public: relation_less_or_equal, relation_greater, relation_greater_or_equal, - + // answered by class symbol symbol, - + // answered by class lst list, - + // answered by class exprseq exprseq, - + // answered by classes numeric, symbol, add, mul, power polynomial, integer_polynomial, @@ -88,17 +88,18 @@ public: rational_polynomial, crational_polynomial, rational_function, - + algebraic, + // answered by class indexed indexed, // class can carry indices has_indices, // object has at least one index - + // answered by class idx idx, - + // answered by class coloridx coloridx, - + // answered by class lorentzidx lorentzidx }; diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp index f76aa9d7..c5d056c2 100644 --- a/ginac/matrix.cpp +++ b/ginac/matrix.cpp @@ -460,24 +460,28 @@ ex matrix::determinant(void) const bool normal_flag = false; unsigned sparse_count = 0; // count non-zero elements for (exvector::const_iterator r=m.begin(); r!=m.end(); ++r) { - if (!(*r).is_zero()) { + if (!(*r).is_zero()) ++sparse_count; - } - if (!(*r).info(info_flags::numeric)) { + if (!(*r).info(info_flags::numeric)) numeric_flag = false; - } if ((*r).info(info_flags::rational_function) && - !(*r).info(info_flags::crational_polynomial)) { + !(*r).info(info_flags::crational_polynomial)) normal_flag = true; - } } if (numeric_flag) // purely numeric matrix return determinant_numeric(); - // Does anybody really know when a matrix is sparse? - if (4*sparse_count0) - sign = -sign; - if (r1>0) { - divisor = tmp.m[(r1-1)*col+(r1-1)].expand(); - // delete the elements we don't need anymore: - for (unsigned c=0; c1, where it can be shown by means of the + // Sylvester determinant that this really divides m[k+1](r,c). + // + // We also allow rational functions where the original prove still holds. + // However, we must care for numerator and denominator separately and + // "manually" work in the integral domains because of subtle cancellations + // (see below). This blows up the bookkeeping a bit and the formula has + // to be modified to expand like this (N{x} stands for numerator of x, + // D{x} for denominator of x): + // N{m[k+1](r,c)} = N{m[k](k,k)}*N{m[k](r,c)}*D{m[k](r,k)}*D{m[k](k,c)} + // -N{m[k](r,k)}*N{m[k](k,c)}*D{m[k](k,k)}*D{m[k](r,c)} + // D{m[k+1](r,c)} = D{m[k](k,k)}*D{m[k](r,c)}*D{m[k](r,k)}*D{m[k](k,c)} + // where for k>1 we now divide N{m[k+1](r,c)} by + // N{m[k-1](k-1,k-1)} + // and D{m[k+1](r,c)} by + // D{m[k-1](k-1,k-1)}. - // first normal all elements: - for (exvector::iterator i=m.begin(); i!=m.end(); ++i) - (*i) = (*i).normal(); + GINAC_ASSERT(det || row==col); + ensure_if_modifiable(); + if (rows()==1) + return 1; - // FIXME: this is unfinished, once matrix::determinant_bareiss is - // bulletproof, some code ought to be copy from there to here. int sign = 1; - ex divisor = 1; - ex dividend; - lst srl; // symbol replacement list for .to_rational() + ex divisor_n = 1; + ex divisor_d = 1; + ex dividend_n; + ex dividend_d; + + // We populate temporary matrices to subsequently operate on. There is + // one holding numerators and another holding denominators of entries. + // This is a must since the evaluator (or even earlier mul's constructor) + // might cancel some trivial element which causes divide() to fail. The + // elements are normalized first (yes, even though this algorithm doesn't + // need GCDs) since the elements of *this might be unnormalized, which + // makes things more complicated than they need to be. + matrix tmp_n(*this); + matrix tmp_d(row,col); // for denominators, if needed + lst srl; // symbol replacement list + exvector::iterator it = m.begin(); + exvector::iterator tmp_n_it = tmp_n.m.begin(); + exvector::iterator tmp_d_it = tmp_d.m.begin(); + for (; it!= m.end(); ++it, ++tmp_n_it, ++tmp_d_it) { + (*tmp_n_it) = (*it).normal().to_rational(srl); + (*tmp_d_it) = (*tmp_n_it).denom(); + (*tmp_n_it) = (*tmp_n_it).numer(); + } for (unsigned r1=0; r10) +// cout << "==<" << r1 << ">" << string(60,'=') << endl; + int indx = tmp_n.pivot(r1); + if (det && indx==-1) + return 0; // FIXME: what to do if det is false? + if (indx>0) { sign = -sign; - if (r1>0) - divisor = this->m[(r1-1)*col+(r1-1)].expand(); + // rows r1 and indx were swapped, so pivot matrix tmp_d: + for (unsigned c=0; c0) { + divisor_n = tmp_n.m[(r1-1)*col+(r1-1)].expand(); + divisor_d = tmp_d.m[(r1-1)*col+(r1-1)].expand(); + // save space by deleting no longer needed elements: + if (det) { + for (unsigned c=0; cm[r1*col+r1]*this->m[r2*col+c] - -this->m[r2*col+r1]*this->m[r1*col+c]).expand(); -#ifdef DO_GINAC_ASSERT - GINAC_ASSERT(divide(dividend.to_rational(srl), - divisor.to_rational(srl), - this->m[r2*col+c])); -#else - divide(dividend.to_rational(srl), - divisor.to_rational(srl), - this->m[r2*col+c]); -#endif // DO_GINAC_ASSERT - this->m[r2*col+c] = this->m[r2*col+c].subs(srl); + dividend_n = (tmp_n.m[r1*col+r1]*tmp_n.m[r2*col+c]* + tmp_d.m[r2*col+r1]*tmp_d.m[r1*col+c] + -tmp_n.m[r2*col+r1]*tmp_n.m[r1*col+c]* + tmp_d.m[r1*col+r1]*tmp_d.m[r2*col+c]).expand(); + dividend_d = (tmp_d.m[r2*col+r1]*tmp_d.m[r1*col+c]* + tmp_d.m[r1*col+r1]*tmp_d.m[r2*col+c]).expand(); +// cout << "Element " << r2 << ',' << c << endl; +// cout << "dividend_n==" << dividend_n << endl; +// cout << "dividend_d==" << dividend_d << endl; +// cout << " divisor_n==" << divisor_n << endl; +// cout << " divisor_d==" << divisor_d << endl; +// cout << string(20,'-') << endl; + bool check = divide(dividend_n, divisor_n, + tmp_n.m[r2*col+c],true); + check &= divide(dividend_d, divisor_d, + tmp_d.m[r2*col+c],true); + GINAC_ASSERT(check); } + // fill up left hand side. for (unsigned c=0; c<=r1; ++c) - this->m[r2*col+c] = _ex0(); + tmp_n.m[r2*col+c] = _ex0(); } +// cout << tmp_n << endl; +// cout << tmp_d << endl; } + // repopulate *this matrix: + it = m.begin(); + tmp_n_it = tmp_n.m.begin(); + tmp_d_it = tmp_d.m.begin(); + for (; it!= m.end(); ++it, ++tmp_n_it, ++tmp_d_it) + (*it) = ((*tmp_n_it)/(*tmp_d_it)).subs(srl); return sign; } diff --git a/ginac/matrix.h b/ginac/matrix.h index a36276f4..9742f398 100644 --- a/ginac/matrix.h +++ b/ginac/matrix.h @@ -97,10 +97,9 @@ public: protected: ex determinant_numeric(void) const; ex determinant_minor(void) const; - ex determinant_bareiss(bool normalize=true) const; int gauss_elimination(void); - int fraction_free_elimination(void); int division_free_elimination(void); + int fraction_free_elimination(bool det = false); int pivot(unsigned ro, bool symbolic=true); private: // FIXME: these should be obsoleted void ffe_swap(unsigned r1, unsigned c1, unsigned r2 ,unsigned c2); diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 860fb10a..e6d4a617 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -280,21 +280,28 @@ void mul::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) cons bool mul::info(unsigned inf) const { - // TODO: optimize - if (inf==info_flags::polynomial || - inf==info_flags::integer_polynomial || - inf==info_flags::cinteger_polynomial || - inf==info_flags::rational_polynomial || - inf==info_flags::crational_polynomial || - inf==info_flags::rational_function) { - for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) { - if (!(recombine_pair_to_ex(*it).info(inf))) - return false; + switch (inf) { + case info_flags::polynomial: + case info_flags::integer_polynomial: + case info_flags::cinteger_polynomial: + case info_flags::rational_polynomial: + case info_flags::crational_polynomial: + case info_flags::rational_function: { + for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) { + if (!(recombine_pair_to_ex(*i).info(inf))) + return false; + } + return overall_coeff.info(inf); + } + case info_flags::algebraic: { + for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) { + if ((recombine_pair_to_ex(*i).info(inf))) + return true; + } + return false; } - return overall_coeff.info(inf); - } else { - return inherited::info(inf); } + return inherited::info(inf); } typedef vector intvector; @@ -419,7 +426,7 @@ ex mul::eval(int level) const ex_to_numeric(addref.overall_coeff). mul_dyn(ex_to_numeric(overall_coeff)))) ->setflag(status_flags::dynallocated | - status_flags::evaluated ); + status_flags::evaluated); } return this->hold(); } diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 4e814361..f93ecc86 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -516,7 +516,8 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args) return true; } #endif - if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial))) + if (check_args && (!a.info(info_flags::rational_polynomial) || + !b.info(info_flags::rational_polynomial))) throw(std::invalid_argument("divide: arguments must be polynomials over the rationals")); // Find first symbol diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp index 49e65be4..3652f4d4 100644 --- a/ginac/numeric.cpp +++ b/ginac/numeric.cpp @@ -221,12 +221,12 @@ numeric::numeric(double d) : basic(TINFO_numeric) numeric::numeric(const char *s) : basic(TINFO_numeric) -{ // MISSING: treatment of complex and ints and rationals. +{ // MISSING: treatment of complex numbers debugmsg("numeric constructor from string",LOGLEVEL_CONSTRUCT); if (strchr(s, '.')) value = new cl_LF(s); else - value = new cl_I(s); + value = new cl_R(s); calchash(); setflag(status_flags::evaluated| status_flags::hash_calculated); @@ -489,42 +489,44 @@ void numeric::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) bool numeric::info(unsigned inf) const { switch (inf) { - case info_flags::numeric: - case info_flags::polynomial: - case info_flags::rational_function: - return true; - case info_flags::real: - return is_real(); - case info_flags::rational: - case info_flags::rational_polynomial: - return is_rational(); - case info_flags::crational: - case info_flags::crational_polynomial: - return is_crational(); - case info_flags::integer: - case info_flags::integer_polynomial: - return is_integer(); - case info_flags::cinteger: - case info_flags::cinteger_polynomial: - return is_cinteger(); - case info_flags::positive: - return is_positive(); - case info_flags::negative: - return is_negative(); - case info_flags::nonnegative: - return !is_negative(); - case info_flags::posint: - return is_pos_integer(); - case info_flags::negint: - return is_integer() && is_negative(); - case info_flags::nonnegint: - return is_nonneg_integer(); - case info_flags::even: - return is_even(); - case info_flags::odd: - return is_odd(); - case info_flags::prime: - return is_prime(); + case info_flags::numeric: + case info_flags::polynomial: + case info_flags::rational_function: + return true; + case info_flags::real: + return is_real(); + case info_flags::rational: + case info_flags::rational_polynomial: + return is_rational(); + case info_flags::crational: + case info_flags::crational_polynomial: + return is_crational(); + case info_flags::integer: + case info_flags::integer_polynomial: + return is_integer(); + case info_flags::cinteger: + case info_flags::cinteger_polynomial: + return is_cinteger(); + case info_flags::positive: + return is_positive(); + case info_flags::negative: + return is_negative(); + case info_flags::nonnegative: + return !is_negative(); + case info_flags::posint: + return is_pos_integer(); + case info_flags::negint: + return is_integer() && is_negative(); + case info_flags::nonnegint: + return is_nonneg_integer(); + case info_flags::even: + return is_even(); + case info_flags::odd: + return is_odd(); + case info_flags::prime: + return is_prime(); + case info_flags::algebraic: + return !is_real(); } return false; } diff --git a/ginac/power.cpp b/ginac/power.cpp index fb4edd56..8ec1945f 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -250,17 +250,20 @@ void power::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) co bool power::info(unsigned inf) const { - if (inf==info_flags::polynomial || - inf==info_flags::integer_polynomial || - inf==info_flags::cinteger_polynomial || - inf==info_flags::rational_polynomial || - inf==info_flags::crational_polynomial) { - return exponent.info(info_flags::nonnegint); - } else if (inf==info_flags::rational_function) { - return exponent.info(info_flags::integer); - } else { - return inherited::info(inf); + switch (inf) { + case info_flags::polynomial: + case info_flags::integer_polynomial: + case info_flags::cinteger_polynomial: + case info_flags::rational_polynomial: + case info_flags::crational_polynomial: + return exponent.info(info_flags::nonnegint); + case info_flags::rational_function: + return exponent.info(info_flags::integer); + case info_flags::algebraic: + return (!exponent.info(info_flags::integer) || + basis.info(inf)); } + return inherited::info(inf); } unsigned power::nops() const