Fix mul::conjugate().
authorRichard Kreckel <kreckel@ginac.de>
Sun, 15 May 2011 16:40:44 +0000 (18:40 +0200)
committerRichard Kreckel <kreckel@ginac.de>
Sun, 15 May 2011 16:40:44 +0000 (18:40 +0200)
Overwrite conjugate() for mul objects because the base class's conjugate
function is incorrect at branch cuts when exponents are fractional (roots).
This has been reported as Sage #10964.

While at it, fix some funny indentation and remove unneeded tests.

check/exam_paranoia.cpp
ginac/expairseq.cpp
ginac/expairseq.h
ginac/inifcns_trans.cpp
ginac/mul.cpp
ginac/mul.h
ginac/pseries.cpp

index 292df81..5fdfb35 100644 (file)
@@ -480,6 +480,20 @@ static unsigned exam_paranoia18()
        return 0;
 }
 
        return 0;
 }
 
+// Bug in mul::conjugate when factors are evaluated at branch cuts, reported as
+// Sage bug #10964.
+static unsigned exam_paranoia19()
+{
+       symbol a("a");
+       ex e = conjugate(a*sqrt(ex(-2))*sqrt(ex(-3)));
+       ex c = a*conjugate(sqrt(ex(-2)))*conjugate(sqrt(ex(-3)));
+       if (!subs(e-c, a==42).is_zero()) {
+               clog << "subs(a*conjugate(sqrt(-2))*conjugate(sqrt(-3))-conjugate(a*sqrt(-2)*sqrt(-3)),a==42) failed to evaluate to 0\n";
+               return 1;
+       }
+       return 0;
+}
+
 unsigned exam_paranoia()
 {
        unsigned result = 0;
 unsigned exam_paranoia()
 {
        unsigned result = 0;
@@ -504,6 +518,7 @@ unsigned exam_paranoia()
        result += exam_paranoia16();  cout << '.' << flush;
        result += exam_paranoia17();  cout << '.' << flush;
        result += exam_paranoia18();  cout << '.' << flush;
        result += exam_paranoia16();  cout << '.' << flush;
        result += exam_paranoia17();  cout << '.' << flush;
        result += exam_paranoia18();  cout << '.' << flush;
+       result += exam_paranoia19();  cout << '.' << flush;
        
        return result;
 }
        
        return result;
 }
index 24991a7..90ec763 100644 (file)
@@ -371,9 +371,7 @@ ex expairseq::conjugate() const
                return *this;
        }
        ex result = thisexpairseq(newepv ? *newepv : seq, x);
                return *this;
        }
        ex result = thisexpairseq(newepv ? *newepv : seq, x);
-       if (newepv) {
-               delete newepv;
-       }
+       delete newepv;
        return result;
 }
 
        return result;
 }
 
index 3fba416..7e8d551 100644 (file)
@@ -108,9 +108,9 @@ protected:
                               unsigned upper_precedence) const;
        virtual expair split_ex_to_pair(const ex & e) const;
        virtual expair combine_ex_with_coeff_to_pair(const ex & e,
                               unsigned upper_precedence) const;
        virtual expair split_ex_to_pair(const ex & e) const;
        virtual expair combine_ex_with_coeff_to_pair(const ex & e,
-                                                                                                const ex & c) const;
+                                                    const ex & c) const;
        virtual expair combine_pair_with_coeff_to_pair(const expair & p,
        virtual expair combine_pair_with_coeff_to_pair(const expair & p,
-                                                                                                  const ex & c) const;
+                                                      const ex & c) const;
        virtual ex recombine_pair_to_ex(const expair & p) const;
        virtual bool expair_needs_further_processing(epp it);
        virtual ex default_overall_coeff() const;
        virtual ex recombine_pair_to_ex(const expair & p) const;
        virtual bool expair_needs_further_processing(epp it);
        virtual ex default_overall_coeff() const;
index da8b3b0..10a3675 100644 (file)
@@ -1249,7 +1249,7 @@ static ex tanh_imag_part(const ex & x)
 
 static ex tanh_conjugate(const ex & x)
 {
 
 static ex tanh_conjugate(const ex & x)
 {
-       // conjugate(tan(x))==tan(conjugate(x))
+       // conjugate(tanh(x))==tanh(conjugate(x))
        return tanh(x.conjugate());
 }
 
        return tanh(x.conjugate());
 }
 
index b8ba32e..6810cb8 100644 (file)
@@ -831,6 +831,38 @@ retry1:
        return ((*this)/divide_by)*multiply_by;
 }
 
        return ((*this)/divide_by)*multiply_by;
 }
 
+ex mul::conjugate() const
+{
+       // The base class' method is wrong here because we have to be careful at
+       // branch cuts. power::conjugate takes care of that already, so use it.
+       epvector *newepv = 0;
+       for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
+               if (newepv) {
+                       newepv->push_back(split_ex_to_pair(recombine_pair_to_ex(*i).conjugate()));
+                       continue;
+               }
+               ex x = recombine_pair_to_ex(*i);
+               ex c = x.conjugate();
+               if (c.is_equal(x)) {
+                       continue;
+               }
+               newepv = new epvector;
+               newepv->reserve(seq.size());
+               for (epvector::const_iterator j=seq.begin(); j!=i; ++j) {
+                       newepv->push_back(*j);
+               }
+               newepv->push_back(split_ex_to_pair(c));
+       }
+       ex x = overall_coeff.conjugate();
+       if (!newepv && are_ex_trivially_equal(x, overall_coeff)) {
+               return *this;
+       }
+       ex result = thisexpairseq(newepv ? *newepv : seq, x);
+       delete newepv;
+       return result;
+}
+
+
 // protected
 
 /** Implementation of ex::diff() for a product.  It applies the product rule.
 // protected
 
 /** Implementation of ex::diff() for a product.  It applies the product rule.
index 0222263..65f59bd 100644 (file)
@@ -64,6 +64,7 @@ public:
        ex smod(const numeric &xi) const;
        numeric max_coefficient() const;
        exvector get_free_indices() const;
        ex smod(const numeric &xi) const;
        numeric max_coefficient() const;
        exvector get_free_indices() const;
+       ex conjugate() const;
 protected:
        ex derivative(const symbol & s) const;
        ex eval_ncmul(const exvector & v) const;
 protected:
        ex derivative(const symbol & s) const;
        ex eval_ncmul(const exvector & v) const;
index f95bbf6..4376817 100644 (file)
@@ -424,14 +424,12 @@ ex pseries::conjugate() const
        epvector * newseq = conjugateepvector(seq);
        ex newpoint = point.conjugate();
 
        epvector * newseq = conjugateepvector(seq);
        ex newpoint = point.conjugate();
 
-       if (!newseq     && are_ex_trivially_equal(point, newpoint)) {
+       if (!newseq && are_ex_trivially_equal(point, newpoint)) {
                return *this;
        }
 
        ex result = (new pseries(var==newpoint, newseq ? *newseq : seq))->setflag(status_flags::dynallocated);
                return *this;
        }
 
        ex result = (new pseries(var==newpoint, newseq ? *newseq : seq))->setflag(status_flags::dynallocated);
-       if (newseq) {
-               delete newseq;
-       }
+       delete newseq;
        return result;
 }
 
        return result;
 }