]> www.ginac.de Git - ginac.git/blobdiff - ginac/power.cpp
synced to 1.1 (expand() problem)
[ginac.git] / ginac / power.cpp
index 0d24a04f950f312c54934ae5b8d2a0e994e7e117..a45bf6afaf02aea954ebd6f2d6c0fe44476b3465 100644 (file)
 #include "indexed.h"
 #include "symbol.h"
 #include "lst.h"
-#include "print.h"
 #include "archive.h"
 #include "utils.h"
 
 namespace GiNaC {
 
-GINAC_IMPLEMENT_REGISTERED_CLASS(power, basic)
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(power, basic,
+  print_func<print_dflt>(&power::do_print_dflt).
+  print_func<print_latex>(&power::do_print_latex).
+  print_func<print_csrc>(&power::do_print_csrc).
+  print_func<print_python>(&power::do_print_python).
+  print_func<print_python_repr>(&power::do_print_python_repr))
 
 typedef std::vector<int> intvector;
 
@@ -85,11 +89,58 @@ DEFAULT_UNARCHIVE(power)
 
 // public
 
+void power::print_power(const print_context & c, const char *powersymbol, const char *openbrace, const char *closebrace, unsigned level) const
+{
+       // Ordinary output of powers using '^' or '**'
+       if (precedence() <= level)
+               c.s << openbrace << '(';
+       basis.print(c, precedence());
+       c.s << powersymbol;
+       c.s << openbrace;
+       exponent.print(c, precedence());
+       c.s << closebrace;
+       if (precedence() <= level)
+               c.s << ')' << closebrace;
+}
+
+void power::do_print_dflt(const print_dflt & c, unsigned level) const
+{
+       if (exponent.is_equal(_ex1_2)) {
+
+               // Square roots are printed in a special way
+               c.s << "sqrt(";
+               basis.print(c);
+               c.s << ')';
+
+       } else
+               print_power(c, "^", "", "", level);
+}
+
+void power::do_print_latex(const print_latex & c, unsigned level) const
+{
+       if (is_exactly_a<numeric>(exponent) && ex_to<numeric>(exponent).is_negative()) {
+
+               // Powers with negative numeric exponents are printed as fractions
+               c.s << "\\frac{1}{";
+               power(basis, -exponent).eval().print(c);
+               c.s << '}';
+
+       } else if (exponent.is_equal(_ex1_2)) {
+
+               // Square roots are printed in a special way
+               c.s << "\\sqrt{";
+               basis.print(c);
+               c.s << '}';
+
+       } else
+               print_power(c, "^", "{", "}", level);
+}
+
 static void print_sym_pow(const print_context & c, const symbol &x, int exp)
 {
        // Optimal output of integer powers of symbols to aid compiler CSE.
        // C.f. ISO/IEC 14882:1998, section 1.9 [intro execution], paragraph 15
-       // to learn why such a parenthisation is really necessary.
+       // to learn why such a parenthesation is really necessary.
        if (exp == 1) {
                x.print(c);
        } else if (exp == 2) {
@@ -109,96 +160,58 @@ static void print_sym_pow(const print_context & c, const symbol &x, int exp)
        }
 }
 
-void power::print(const print_context & c, unsigned level) const
+void power::do_print_csrc(const print_csrc & c, unsigned level) const
 {
-       if (is_a<print_tree>(c)) {
-
-               inherited::print(c, level);
-
-       } else if (is_a<print_csrc>(c)) {
-
-               // Integer powers of symbols are printed in a special, optimized way
-               if (exponent.info(info_flags::integer)
-                && (is_a<symbol>(basis) || is_a<constant>(basis))) {
-                       int exp = ex_to<numeric>(exponent).to_int();
-                       if (exp > 0)
-                               c.s << '(';
-                       else {
-                               exp = -exp;
-                               if (is_a<print_csrc_cl_N>(c))
-                                       c.s << "recip(";
-                               else
-                                       c.s << "1.0/(";
-                       }
-                       print_sym_pow(c, ex_to<symbol>(basis), exp);
-                       c.s << ')';
-
-               // <expr>^-1 is printed as "1.0/<expr>" or with the recip() function of CLN
-               } else if (exponent.is_equal(_ex_1)) {
+       // Integer powers of symbols are printed in a special, optimized way
+       if (exponent.info(info_flags::integer)
+        && (is_a<symbol>(basis) || is_a<constant>(basis))) {
+               int exp = ex_to<numeric>(exponent).to_int();
+               if (exp > 0)
+                       c.s << '(';
+               else {
+                       exp = -exp;
                        if (is_a<print_csrc_cl_N>(c))
                                c.s << "recip(";
                        else
                                c.s << "1.0/(";
-                       basis.print(c);
-                       c.s << ')';
-
-               // Otherwise, use the pow() or expt() (CLN) functions
-               } else {
-                       if (is_a<print_csrc_cl_N>(c))
-                               c.s << "expt(";
-                       else
-                               c.s << "pow(";
-                       basis.print(c);
-                       c.s << ',';
-                       exponent.print(c);
-                       c.s << ')';
                }
+               print_sym_pow(c, ex_to<symbol>(basis), exp);
+               c.s << ')';
 
-       } else if (is_a<print_python_repr>(c)) {
+       // <expr>^-1 is printed as "1.0/<expr>" or with the recip() function of CLN
+       } else if (exponent.is_equal(_ex_1)) {
+               if (is_a<print_csrc_cl_N>(c))
+                       c.s << "recip(";
+               else
+                       c.s << "1.0/(";
+               basis.print(c);
+               c.s << ')';
 
-               c.s << class_name() << '(';
+       // Otherwise, use the pow() or expt() (CLN) functions
+       } else {
+               if (is_a<print_csrc_cl_N>(c))
+                       c.s << "expt(";
+               else
+                       c.s << "pow(";
                basis.print(c);
                c.s << ',';
                exponent.print(c);
                c.s << ')';
+       }
+}
 
-       } else {
-
-               bool is_tex = is_a<print_latex>(c);
-
-               if (is_tex && is_exactly_a<numeric>(exponent) && ex_to<numeric>(exponent).is_negative()) {
-
-                       // Powers with negative numeric exponents are printed as fractions in TeX
-                       c.s << "\\frac{1}{";
-                       power(basis, -exponent).eval().print(c);
-                       c.s << "}";
-
-               } else if (exponent.is_equal(_ex1_2)) {
-
-                       // Square roots are printed in a special way
-                       c.s << (is_tex ? "\\sqrt{" : "sqrt(");
-                       basis.print(c);
-                       c.s << (is_tex ? '}' : ')');
-
-               } else {
+void power::do_print_python(const print_python & c, unsigned level) const
+{
+       print_power(c, "**", "", "", level);
+}
 
-                       // Ordinary output of powers using '^' or '**'
-                       if (precedence() <= level)
-                               c.s << (is_tex ? "{(" : "(");
-                       basis.print(c, precedence());
-                       if (is_a<print_python>(c))
-                               c.s << "**";
-                       else
-                               c.s << '^';
-                       if (is_tex)
-                               c.s << '{';
-                       exponent.print(c, precedence());
-                       if (is_tex)
-                               c.s << '}';
-                       if (precedence() <= level)
-                               c.s << (is_tex ? ")}" : ")");
-               }
-       }
+void power::do_print_python_repr(const print_python_repr & c, unsigned level) const
+{
+       c.s << class_name() << '(';
+       basis.print(c);
+       c.s << ',';
+       exponent.print(c);
+       c.s << ')';
 }
 
 bool power::info(unsigned inf) const
@@ -518,27 +531,26 @@ ex power::evalm() const
 // from mul.cpp
 extern bool tryfactsubs(const ex &, const ex &, int &, lst &);
 
-ex power::subs(const lst & ls, const lst & lr, unsigned options) const
+ex power::subs(const exmap & m, unsigned options) const
 {      
-       const ex &subsed_basis = basis.subs(ls, lr, options);
-       const ex &subsed_exponent = exponent.subs(ls, lr, options);
+       const ex &subsed_basis = basis.subs(m, options);
+       const ex &subsed_exponent = exponent.subs(m, options);
 
        if (!are_ex_trivially_equal(basis, subsed_basis)
         || !are_ex_trivially_equal(exponent, subsed_exponent)) 
-               return power(subsed_basis, subsed_exponent).subs_one_level(ls, lr, options);
+               return power(subsed_basis, subsed_exponent).subs_one_level(m, options);
 
-       if (!(options & subs_options::subs_algebraic))
-               return subs_one_level(ls, lr, options);
+       if (!(options & subs_options::algebraic))
+               return subs_one_level(m, options);
 
-       lst::const_iterator its, itr;
-       for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
+       for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) {
                int nummatches = std::numeric_limits<int>::max();
                lst repls;
-               if (tryfactsubs(*this, *its, nummatches, repls))
-                       return (ex_to<basic>((*this) * power(itr->subs(ex(repls), subs_options::subs_no_pattern) / its->subs(ex(repls), subs_options::subs_no_pattern), nummatches))).subs_one_level(ls, lr, options);
+               if (tryfactsubs(*this, it->first, nummatches, repls))
+                       return (ex_to<basic>((*this) * power(it->second.subs(ex(repls), subs_options::no_pattern) / it->first.subs(ex(repls), subs_options::no_pattern), nummatches))).subs_one_level(m, options);
        }
 
-       return subs_one_level(ls, lr, options);
+       return subs_one_level(m, options);
 }
 
 ex power::eval_ncmul(const exvector & v) const
@@ -827,6 +839,8 @@ ex power::expand_mul(const mul & m, const numeric & n) const
 
        epvector distrseq;
        distrseq.reserve(m.seq.size());
+       bool need_reexpand = false;
+
        epvector::const_iterator last = m.seq.end();
        epvector::const_iterator cit = m.seq.begin();
        while (cit!=last) {
@@ -835,11 +849,22 @@ ex power::expand_mul(const mul & m, const numeric & n) const
                } else {
                        // it is safe not to call mul::combine_pair_with_coeff_to_pair()
                        // since n is an integer
-                       distrseq.push_back(expair(cit->rest, ex_to<numeric>(cit->coeff).mul(n)));
+                       numeric new_coeff = ex_to<numeric>(cit->coeff).mul(n);
+                       if (is_exactly_a<add>(cit->rest) && new_coeff.is_pos_integer()) {
+                               // this happens when e.g. (a+b)^(1/2) gets squared and
+                               // the resulting product needs to be reexpanded
+                               need_reexpand = true;
+                       }
+                       distrseq.push_back(expair(cit->rest, new_coeff));
                }
                ++cit;
        }
-       return (new mul(distrseq, ex_to<numeric>(m.overall_coeff).power_dyn(n)))->setflag(status_flags::dynallocated);
+
+       const mul & result = static_cast<const mul &>((new mul(distrseq, ex_to<numeric>(m.overall_coeff).power_dyn(n)))->setflag(status_flags::dynallocated));
+       if (need_reexpand)
+               return ex(result).expand();
+       else
+               return result.setflag(status_flags::expanded);
 }
 
 } // namespace GiNaC