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.

diff --git a/check/exam_paranoia.cpp b/check/exam_paranoia.cpp
index 61bd85debb427dd380b9c1fa3145a0292f140e44..9e6fb015f6b84e05be483d9dde9155e6df62a8e2 100644 (file)
--- 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();
 }
 
        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()
 {
 
 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_paranoia15();  cout << '.' << flush;
        result += exam_paranoia16();  cout << '.' << flush;
        result += exam_paranoia17();  cout << '.' << flush;

+       result += exam_paranoia18();  cout << '.' << flush;

        
        return result;
 }
        
        return result;
 }

diff --git a/ginac/add.cpp b/ginac/add.cpp
index 6548366ee16f1facb3dab2c1751444730171f236..241516e07283b7e49a3a7fd7bf47490bebb1c9f9 100644 (file)
--- a/ginac/add.cpp
+++ b/ginac/add.cpp
@@ -28,6 +28,7 @@
 #include "utils.h"
 #include "clifford.h"
 #include "ncmul.h"
 #include "utils.h"
 #include "clifford.h"
 #include "ncmul.h"

+#include "compiler.h"

 
 #include <iostream>
 #include <limits>
 
 #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));
        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
                ++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"));
        }
        } 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();
 }
 
        return this->hold();
 }