From 89d5356b4aa33cb4481575f9453f36c3404b015b Mon Sep 17 00:00:00 2001
From: Richard Kreckel <kreckel@ginac.de>
Date: Thu, 23 Sep 2010 00:40:38 +0200
Subject: [PATCH] Make sure add::eval() collects all numeric terms.

Apparently, add::eval() assumed that none of the elements of its epvector
has a numeric rest. However, nothing guarantees that -- in particular
evalchildren() doesn't (and actually cannot) do so. Since there are many
places where a new add is constructed directly from an epvector, enforcing
this doesn't make sense either. One example where it did fail was found by
Burgin Erocal: real_part(1+2*(sqrt(2)+1)*(sqrt(2)-1)) returned 1+2*1, not 3.

Thanks to Burcin Erocal for reporting this bug.
(cherry picked from commit e08cda1854bdb82f6706ec269233577690ae00e4)
---
 check/exam_paranoia.cpp | 13 +++++++++++++
 ginac/add.cpp           | 31 ++++++++++++++++++++++++++++---
 2 files changed, 41 insertions(+), 3 deletions(-)

diff --git a/check/exam_paranoia.cpp b/check/exam_paranoia.cpp
index 61bd85de..9e6fb015 100644
--- a/check/exam_paranoia.cpp
+++ b/check/exam_paranoia.cpp
@@ -467,6 +467,18 @@ static unsigned exam_paranoia17()
 	return test_cycl.get_free_indices().size();
 }
 
+// Bug in add::eval() could result in numeric terms not being collected into
+// the overall coefficient. Fixed on Sep 22, 2010
+static unsigned exam_paranoia18()
+{
+	ex sqrt2 = sqrt(ex(2));
+	ex e = 1+2*(sqrt2+1)*(sqrt2-1);
+	if ( e.real_part() != 3 ) {
+		clog << "real_part(1+2*(sqrt(2)+1)*(sqrt(2)-1)) failed to evaluate to 3\n";
+		return 1;
+	}
+	return 0;
+}
 
 unsigned exam_paranoia()
 {
@@ -491,6 +503,7 @@ unsigned exam_paranoia()
 	result += exam_paranoia15();  cout << '.' << flush;
 	result += exam_paranoia16();  cout << '.' << flush;
 	result += exam_paranoia17();  cout << '.' << flush;
+	result += exam_paranoia18();  cout << '.' << flush;
 	
 	return result;
 }
diff --git a/ginac/add.cpp b/ginac/add.cpp
index 6548366e..241516e0 100644
--- a/ginac/add.cpp
+++ b/ginac/add.cpp
@@ -28,6 +28,7 @@
 #include "utils.h"
 #include "clifford.h"
 #include "ncmul.h"
+#include "compiler.h"
 
 #include <iostream>
 #include <limits>
@@ -343,9 +344,6 @@ ex add::eval(int level) const
 	epvector::const_iterator i = seq.begin(), end = seq.end();
 	while (i != end) {
 		GINAC_ASSERT(!is_exactly_a<add>(i->rest));
-		if (is_exactly_a<numeric>(i->rest))
-			dbgprint();
-		GINAC_ASSERT(!is_exactly_a<numeric>(i->rest));
 		++i;
 	}
 #endif // def DO_GINAC_ASSERT
@@ -366,6 +364,33 @@ ex add::eval(int level) const
 	} else if (!overall_coeff.is_zero() && seq[0].rest.return_type() != return_types::commutative) {
 		throw (std::logic_error("add::eval(): sum of non-commutative objects has non-zero numeric term"));
 	}
+	
+	// if any terms in the sum still are purely numeric, then they are more
+	// appropriately collected into the overall coefficient
+	epvector::const_iterator last = seq.end();
+	epvector::const_iterator j = seq.begin();
+	int terms_to_collect = 0;
+	while (j != last) {
+		if (unlikely(is_a<numeric>(j->rest)))
+			++terms_to_collect;
+		++j;
+	}
+	if (terms_to_collect) {
+		std::auto_ptr<epvector> s(new epvector);
+		s->reserve(seq_size - terms_to_collect);
+		numeric oc = *_num1_p;
+		j = seq.begin();
+		while (j != last) {
+			if (unlikely(is_a<numeric>(j->rest)))
+				oc = oc.mul(ex_to<numeric>(j->rest)).mul(ex_to<numeric>(j->coeff));
+			else
+				s->push_back(*j);
+			++j;
+		}
+		return (new add(s, ex_to<numeric>(overall_coeff).add_dyn(oc)))
+		        ->setflag(status_flags::dynallocated);
+	}
+	
 	return this->hold();
 }
 
-- 
2.44.0