From f9afb6aca6a971650dff63b11ca8c2ef18506690 Mon Sep 17 00:00:00 2001
From: Richard Kreckel <kreckel@ginac.de>
Date: Sun, 15 May 2011 18:40:44 +0200
Subject: [PATCH] Fix mul::conjugate().

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 | 15 +++++++++++++++
 ginac/expairseq.cpp     |  4 +---
 ginac/expairseq.h       |  4 ++--
 ginac/inifcns_trans.cpp |  2 +-
 ginac/mul.cpp           | 32 ++++++++++++++++++++++++++++++++
 ginac/mul.h             |  1 +
 ginac/pseries.cpp       |  6 ++----
 7 files changed, 54 insertions(+), 10 deletions(-)

diff --git a/check/exam_paranoia.cpp b/check/exam_paranoia.cpp
index 292df819..5fdfb350 100644
--- a/check/exam_paranoia.cpp
+++ b/check/exam_paranoia.cpp
@@ -480,6 +480,20 @@ static unsigned exam_paranoia18()
 	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;
@@ -504,6 +518,7 @@ unsigned exam_paranoia()
 	result += exam_paranoia16();  cout << '.' << flush;
 	result += exam_paranoia17();  cout << '.' << flush;
 	result += exam_paranoia18();  cout << '.' << flush;
+	result += exam_paranoia19();  cout << '.' << flush;
 	
 	return result;
 }
diff --git a/ginac/expairseq.cpp b/ginac/expairseq.cpp
index 24991a77..90ec7635 100644
--- a/ginac/expairseq.cpp
+++ b/ginac/expairseq.cpp
@@ -371,9 +371,7 @@ ex expairseq::conjugate() const
 		return *this;
 	}
 	ex result = thisexpairseq(newepv ? *newepv : seq, x);
-	if (newepv) {
-		delete newepv;
-	}
+	delete newepv;
 	return result;
 }
 
diff --git a/ginac/expairseq.h b/ginac/expairseq.h
index 3fba4161..7e8d5511 100644
--- a/ginac/expairseq.h
+++ b/ginac/expairseq.h
@@ -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,
-												 const ex & c) const;
+	                                             const ex & c) const;
 	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;
diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp
index da8b3b0a..10a36758 100644
--- a/ginac/inifcns_trans.cpp
+++ b/ginac/inifcns_trans.cpp
@@ -1249,7 +1249,7 @@ static ex tanh_imag_part(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());
 }
 
diff --git a/ginac/mul.cpp b/ginac/mul.cpp
index b8ba32ea..6810cb8d 100644
--- a/ginac/mul.cpp
+++ b/ginac/mul.cpp
@@ -831,6 +831,38 @@ retry1:
 	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.
diff --git a/ginac/mul.h b/ginac/mul.h
index 02222637..65f59bde 100644
--- a/ginac/mul.h
+++ b/ginac/mul.h
@@ -64,6 +64,7 @@ public:
 	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;
diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp
index f95bbf6d..43768179 100644
--- a/ginac/pseries.cpp
+++ b/ginac/pseries.cpp
@@ -424,14 +424,12 @@ ex pseries::conjugate() const
 	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);
-	if (newseq) {
-		delete newseq;
-	}
+	delete newseq;
 	return result;
 }
 
-- 
2.44.0