]> www.ginac.de Git - ginac.git/blobdiff - ginac/mul.cpp
synced to 1.1 (expand() problem)
[ginac.git] / ginac / mul.cpp
index 9d2aac1518dfc795e3c440b630884dec8bcad9c6..db68b275ffbbc788a476bad5b9b9f814b30229ff 100644 (file)
 
 namespace GiNaC {
 
-GINAC_IMPLEMENT_REGISTERED_CLASS(mul, expairseq)
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(mul, expairseq,
+  print_func<print_context>(&mul::do_print).
+  print_func<print_latex>(&mul::do_print_latex).
+  print_func<print_csrc>(&mul::do_print_csrc).
+  print_func<print_tree>(&inherited::do_print_tree).
+  print_func<print_python_repr>(&mul::do_print_python_repr))
+
 
 //////////
 // default constructor
@@ -85,13 +91,12 @@ mul::mul(const epvector & v, const ex & oc)
        GINAC_ASSERT(is_canonical());
 }
 
-mul::mul(epvector * vp, const ex & oc)
+mul::mul(std::auto_ptr<epvector> vp, const ex & oc)
 {
        tinfo_key = TINFO_mul;
        GINAC_ASSERT(vp!=0);
        overall_coeff = oc;
        construct_from_epvector(*vp);
-       delete vp;
        GINAC_ASSERT(is_canonical());
 }
 
@@ -118,159 +123,154 @@ DEFAULT_ARCHIVING(mul)
 // functions overriding virtual functions from base classes
 //////////
 
-// public
-void mul::print(const print_context & c, unsigned level) const
+void mul::print_overall_coeff(const print_context & c, const char *mul_sym) const
 {
-       if (is_a<print_tree>(c)) {
+       const numeric &coeff = ex_to<numeric>(overall_coeff);
+       if (coeff.csgn() == -1)
+               c.s << '-';
+       if (!coeff.is_equal(_num1) &&
+               !coeff.is_equal(_num_1)) {
+               if (coeff.is_rational()) {
+                       if (coeff.is_negative())
+                               (-coeff).print(c);
+                       else
+                               coeff.print(c);
+               } else {
+                       if (coeff.csgn() == -1)
+                               (-coeff).print(c, precedence());
+                       else
+                               coeff.print(c, precedence());
+               }
+               c.s << mul_sym;
+       }
+}
 
-               inherited::print(c, level);
+void mul::do_print(const print_context & c, unsigned level) const
+{
+       if (precedence() <= level)
+               c.s << '(';
 
-       } else if (is_a<print_csrc>(c)) {
+       print_overall_coeff(c, "*");
 
-               if (precedence() <= level)
-                       c.s << "(";
+       epvector::const_iterator it = seq.begin(), itend = seq.end();
+       bool first = true;
+       while (it != itend) {
+               if (!first)
+                       c.s << '*';
+               else
+                       first = false;
+               recombine_pair_to_ex(*it).print(c, precedence());
+               ++it;
+       }
 
-               if (!overall_coeff.is_equal(_ex1)) {
-                       overall_coeff.print(c, precedence());
-                       c.s << "*";
-               }
+       if (precedence() <= level)
+               c.s << ')';
+}
 
-               // Print arguments, separated by "*" or "/"
-               epvector::const_iterator it = seq.begin(), itend = seq.end();
-               while (it != itend) {
-
-                       // If the first argument is a negative integer power, it gets printed as "1.0/<expr>"
-                       bool needclosingparenthesis = false;
-                       if (it == seq.begin() && it->coeff.info(info_flags::negint)) {
-                               if (is_a<print_csrc_cl_N>(c)) {
-                                       c.s << "recip(";
-                                       needclosingparenthesis = true;
-                               } else
-                                       c.s << "1.0/";
-                       }
+void mul::do_print_latex(const print_latex & c, unsigned level) const
+{
+       if (precedence() <= level)
+               c.s << "{(";
 
-                       // If the exponent is 1 or -1, it is left out
-                       if (it->coeff.is_equal(_ex1) || it->coeff.is_equal(_ex_1))
-                               it->rest.print(c, precedence());
-                       else if (it->coeff.info(info_flags::negint))
-                               // Outer parens around ex needed for broken GCC parser:
-                               (ex(power(it->rest, -ex_to<numeric>(it->coeff)))).print(c, level);
-                       else
-                               // Outer parens around ex needed for broken GCC parser:
-                               (ex(power(it->rest, ex_to<numeric>(it->coeff)))).print(c, level);
-
-                       if (needclosingparenthesis)
-                               c.s << ")";
-
-                       // Separator is "/" for negative integer powers, "*" otherwise
-                       ++it;
-                       if (it != itend) {
-                               if (it->coeff.info(info_flags::negint))
-                                       c.s << "/";
-                               else
-                                       c.s << "*";
-                       }
-               }
+       print_overall_coeff(c, " ");
 
-               if (precedence() <= level)
-                       c.s << ")";
+       // Separate factors into those with negative numeric exponent
+       // and all others
+       epvector::const_iterator it = seq.begin(), itend = seq.end();
+       exvector neg_powers, others;
+       while (it != itend) {
+               GINAC_ASSERT(is_exactly_a<numeric>(it->coeff));
+               if (ex_to<numeric>(it->coeff).is_negative())
+                       neg_powers.push_back(recombine_pair_to_ex(expair(it->rest, -(it->coeff))));
+               else
+                       others.push_back(recombine_pair_to_ex(*it));
+               ++it;
+       }
 
-       } else if (is_a<print_python_repr>(c)) {
-               c.s << class_name() << '(';
-               op(0).print(c);
-               for (size_t i=1; i<nops(); ++i) {
-                       c.s << ',';
-                       op(i).print(c);
-               }
-               c.s << ')';
-       } else {
+       if (!neg_powers.empty()) {
 
-               if (precedence() <= level) {
-                       if (is_a<print_latex>(c))
-                               c.s << "{(";
-                       else
-                               c.s << "(";
-               }
+               // Factors with negative exponent are printed as a fraction
+               c.s << "\\frac{";
+               mul(others).eval().print(c);
+               c.s << "}{";
+               mul(neg_powers).eval().print(c);
+               c.s << "}";
 
-               // First print the overall numeric coefficient
-               const numeric &coeff = ex_to<numeric>(overall_coeff);
-               if (coeff.csgn() == -1)
-                       c.s << '-';
-               if (!coeff.is_equal(_num1) &&
-                       !coeff.is_equal(_num_1)) {
-                       if (coeff.is_rational()) {
-                               if (coeff.is_negative())
-                                       (-coeff).print(c);
-                               else
-                                       coeff.print(c);
-                       } else {
-                               if (coeff.csgn() == -1)
-                                       (-coeff).print(c, precedence());
-                               else
-                                       coeff.print(c, precedence());
-                       }
-                       if (is_a<print_latex>(c))
-                               c.s << ' ';
-                       else
-                               c.s << '*';
-               }
+       } else {
 
-               // Then proceed with the remaining factors
-               epvector::const_iterator it = seq.begin(), itend = seq.end();
-               if (is_a<print_latex>(c)) {
-
-                       // Separate factors into those with negative numeric exponent
-                       // and all others
-                       exvector neg_powers, others;
-                       while (it != itend) {
-                               GINAC_ASSERT(is_exactly_a<numeric>(it->coeff));
-                               if (ex_to<numeric>(it->coeff).is_negative())
-                                       neg_powers.push_back(recombine_pair_to_ex(expair(it->rest, -(it->coeff))));
-                               else
-                                       others.push_back(recombine_pair_to_ex(*it));
-                               ++it;
-                       }
+               // All other factors are printed in the ordinary way
+               exvector::const_iterator vit = others.begin(), vitend = others.end();
+               while (vit != vitend) {
+                       c.s << ' ';
+                       vit->print(c, precedence());
+                       ++vit;
+               }
+       }
 
-                       if (!neg_powers.empty()) {
+       if (precedence() <= level)
+               c.s << ")}";
+}
 
-                               // Factors with negative exponent are printed as a fraction
-                               c.s << "\\frac{";
-                               mul(others).eval().print(c);
-                               c.s << "}{";
-                               mul(neg_powers).eval().print(c);
-                               c.s << "}";
+void mul::do_print_csrc(const print_csrc & c, unsigned level) const
+{
+       if (precedence() <= level)
+               c.s << "(";
 
-                       } else {
+       if (!overall_coeff.is_equal(_ex1)) {
+               overall_coeff.print(c, precedence());
+               c.s << "*";
+       }
 
-                               // All other factors are printed in the ordinary way
-                               exvector::const_iterator vit = others.begin(), vitend = others.end();
-                               while (vit != vitend) {
-                                       c.s << ' ';
-                                       vit->print(c, precedence());
-                                       ++vit;
-                               }
-                       }
+       // Print arguments, separated by "*" or "/"
+       epvector::const_iterator it = seq.begin(), itend = seq.end();
+       while (it != itend) {
+
+               // If the first argument is a negative integer power, it gets printed as "1.0/<expr>"
+               bool needclosingparenthesis = false;
+               if (it == seq.begin() && it->coeff.info(info_flags::negint)) {
+                       if (is_a<print_csrc_cl_N>(c)) {
+                               c.s << "recip(";
+                               needclosingparenthesis = true;
+                       } else
+                               c.s << "1.0/";
+               }
 
-               } else {
+               // If the exponent is 1 or -1, it is left out
+               if (it->coeff.is_equal(_ex1) || it->coeff.is_equal(_ex_1))
+                       it->rest.print(c, precedence());
+               else if (it->coeff.info(info_flags::negint))
+                       // Outer parens around ex needed for broken GCC parser:
+                       (ex(power(it->rest, -ex_to<numeric>(it->coeff)))).print(c, level);
+               else
+                       // Outer parens around ex needed for broken GCC parser:
+                       (ex(power(it->rest, ex_to<numeric>(it->coeff)))).print(c, level);
 
-                       bool first = true;
-                       while (it != itend) {
-                               if (!first)
-                                       c.s << '*';
-                               else
-                                       first = false;
-                               recombine_pair_to_ex(*it).print(c, precedence());
-                               ++it;
-                       }
-               }
+               if (needclosingparenthesis)
+                       c.s << ")";
 
-               if (precedence() <= level) {
-                       if (is_a<print_latex>(c))
-                               c.s << ")}";
+               // Separator is "/" for negative integer powers, "*" otherwise
+               ++it;
+               if (it != itend) {
+                       if (it->coeff.info(info_flags::negint))
+                               c.s << "/";
                        else
-                               c.s << ")";
+                               c.s << "*";
                }
        }
+
+       if (precedence() <= level)
+               c.s << ")";
+}
+
+void mul::do_print_python_repr(const print_python_repr & c, unsigned level) const
+{
+       c.s << class_name() << '(';
+       op(0).print(c);
+       for (size_t i=1; i<nops(); ++i) {
+               c.s << ',';
+               op(i).print(c);
+       }
+       c.s << ')';
 }
 
 bool mul::info(unsigned inf) const
@@ -378,10 +378,10 @@ ex mul::coeff(const ex & s, int n) const
  *  @param level cut-off in recursive evaluation */
 ex mul::eval(int level) const
 {
-       epvector *evaled_seqp = evalchildren(level);
-       if (evaled_seqp) {
+       std::auto_ptr<epvector> evaled_seqp = evalchildren(level);
+       if (evaled_seqp.get()) {
                // do more evaluation later
-               return (new mul(evaled_seqp,overall_coeff))->
+               return (new mul(evaled_seqp, overall_coeff))->
                           setflag(status_flags::dynallocated);
        }
        
@@ -424,7 +424,7 @@ ex mul::eval(int level) const
                   ex_to<numeric>((*seq.begin()).coeff).is_equal(_num1)) {
                // *(+(x,y,...);c) -> +(*(x,c),*(y,c),...) (c numeric(), no powers of +())
                const add & addref = ex_to<add>((*seq.begin()).rest);
-               epvector *distrseq = new epvector();
+               std::auto_ptr<epvector> distrseq(new epvector);
                distrseq->reserve(addref.seq.size());
                epvector::const_iterator i = addref.seq.begin(), end = addref.seq.end();
                while (i != end) {
@@ -447,7 +447,7 @@ ex mul::evalf(int level) const
        if (level==-max_recursion_level)
                throw(std::runtime_error("max recursion level reached"));
        
-       epvector *s = new epvector();
+       std::auto_ptr<epvector> s(new epvector);
        s->reserve(seq.size());
 
        --level;
@@ -470,7 +470,7 @@ ex mul::evalm() const
        // Evaluate children first, look whether there are any matrices at all
        // (there can be either no matrices or one matrix; if there were more
        // than one matrix, it would be a non-commutative product)
-       epvector *s = new epvector;
+       std::auto_ptr<epvector> s(new epvector);
        s->reserve(seq.size());
 
        bool have_matrix = false;
@@ -723,7 +723,7 @@ ex mul::thisexpairseq(const epvector & v, const ex & oc) const
        return (new mul(v, oc))->setflag(status_flags::dynallocated);
 }
 
-ex mul::thisexpairseq(epvector * vp, const ex & oc) const
+ex mul::thisexpairseq(std::auto_ptr<epvector> vp, const ex & oc) const
 {
        return (new mul(vp, oc))->setflag(status_flags::dynallocated);
 }
@@ -827,16 +827,19 @@ bool mul::can_make_flat(const expair & p) const
 ex mul::expand(unsigned options) const
 {
        // First, expand the children
-       epvector * expanded_seqp = expandchildren(options);
-       const epvector & expanded_seq = (expanded_seqp == NULL) ? seq : *expanded_seqp;
+       std::auto_ptr<epvector> expanded_seqp = expandchildren(options);
+       const epvector & expanded_seq = (expanded_seqp.get() ? *expanded_seqp : seq);
 
        // Now, look for all the factors that are sums and multiply each one out
        // with the next one that is found while collecting the factors which are
        // not sums
        int number_of_adds = 0;
        ex last_expanded = _ex1;
+
        epvector non_adds;
        non_adds.reserve(expanded_seq.size());
+       bool non_adds_has_sums = false; // Look for sums or powers of sums in the non_adds (we need this later)
+
        epvector::const_iterator cit = expanded_seq.begin(), last = expanded_seq.end();
        while (cit != last) {
                if (is_exactly_a<add>(cit->rest) &&
@@ -886,7 +889,7 @@ ex mul::expand(unsigned options) const
                                        for (epvector::const_iterator i2=add2begin; i2!=add2end; ++i2) {
                                                // Don't push_back expairs which might have a rest that evaluates to a numeric,
                                                // since that would violate an invariant of expairseq:
-                                               const ex rest = (new mul(i1->rest, i2->rest))->setflag(status_flags::dynallocated);
+                                               const ex rest = ex((new mul(i1->rest, i2->rest))->setflag(status_flags::dynallocated)).expand();
                                                if (is_exactly_a<numeric>(rest))
                                                        oc += ex_to<numeric>(rest).mul(ex_to<numeric>(i1->coeff).mul(ex_to<numeric>(i2->coeff)));
                                                else
@@ -901,25 +904,36 @@ ex mul::expand(unsigned options) const
                                last_expanded = cit->rest;
                        }
                } else {
+                       if (is_exactly_a<add>(cit->rest))
+                               non_adds_has_sums = true;
                        non_adds.push_back(*cit);
                }
                ++cit;
        }
-       if (expanded_seqp)
-               delete expanded_seqp;
-       
+
        // Now the only remaining thing to do is to multiply the factors which
        // were not sums into the "last_expanded" sum
        if (is_exactly_a<add>(last_expanded)) {
                const add & finaladd = ex_to<add>(last_expanded);
-               exvector distrseq;
+
                size_t n = finaladd.nops();
+               exvector distrseq;
                distrseq.reserve(n);
+
                for (size_t i=0; i<n; ++i) {
                        epvector factors = non_adds;
-                       factors.push_back(split_ex_to_pair(finaladd.op(i)));
-                       distrseq.push_back((new mul(factors, overall_coeff))->
-                                           setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0)));
+                       expair new_factor = split_ex_to_pair(finaladd.op(i).expand());
+                       factors.push_back(new_factor);
+
+                       const mul & term = static_cast<const mul &>((new mul(factors, overall_coeff))->setflag(status_flags::dynallocated));
+
+                       // The new term may have sums in it if e.g. a sqrt() of a sum in
+                       // the non_adds meets a sqrt() of a sum in the factor from
+                       // last_expanded. In this case we should re-expand the term.
+                       if (non_adds_has_sums || is_exactly_a<add>(new_factor.rest))
+                               distrseq.push_back(ex(term).expand());
+                       else
+                               distrseq.push_back(term.setflag(options == 0 ? status_flags::expanded : 0));
                }
                return ((new add(distrseq))->
                        setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0)));
@@ -948,7 +962,7 @@ ex mul::expand(unsigned options) const
  *  @see mul::expand()
  *  @return pointer to epvector containing expanded representation or zero
  *  pointer, if sequence is unchanged. */
-epvector * mul::expandchildren(unsigned options) const
+std::auto_ptr<epvector> mul::expandchildren(unsigned options) const
 {
        const epvector::const_iterator last = seq.end();
        epvector::const_iterator cit = seq.begin();
@@ -958,7 +972,7 @@ epvector * mul::expandchildren(unsigned options) const
                if (!are_ex_trivially_equal(factor,expanded_factor)) {
                        
                        // something changed, copy seq, eval and return it
-                       epvector *s = new epvector;
+                       std::auto_ptr<epvector> s(new epvector);
                        s->reserve(seq.size());
                        
                        // copy parts of seq which are known not to have changed
@@ -967,9 +981,11 @@ epvector * mul::expandchildren(unsigned options) const
                                s->push_back(*cit2);
                                ++cit2;
                        }
+
                        // copy first changed element
                        s->push_back(split_ex_to_pair(expanded_factor));
                        ++cit2;
+
                        // copy rest
                        while (cit2!=last) {
                                s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit2).expand(options)));
@@ -980,7 +996,7 @@ epvector * mul::expandchildren(unsigned options) const
                ++cit;
        }
        
-       return 0; // nothing has changed
+       return std::auto_ptr<epvector>(0); // nothing has changed
 }
 
 } // namespace GiNaC