Fix log() error in multiple polylog G.

Tobias Huber reported the following bug:
> G({-0.18784426860496276*I,-0.1579006477353911*I,+0.18784426860496276*I,-0.18784426860496276*I},1);
log_eval(): log(0)

The problem is the sorting operation in G_do_trafo. This operation should
1) put the elements in increasing order of |x|
2) put equal elements next to each other (which then will avoid the log(0) problem).

The current code fails for goal 2) if we have complex numbers of equal absolute
value, one number occuring more than once, in an initial order like in the example
above:
Entries 1,3,4 of the list all have the same absolute value, entries 1 and 4 are
equal. In the sorting operation 1 and 4 should be put next to each other.
Previously, the sorting operation would give 2,1,3,4. What we would like to have
is either 2,1,4,3 or 2,3,1,4. It is clear that the key for the sorting operation
cannot be |x| alone. If |x1|=|x2| we have to use the phase as well.

This patch takes the phase into account in the sorting.

diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp
index 866472f4081a317e23b1cefcf57702b0545f245e..a42c1bd3245895a697ca46108fe206afef175542 100644 (file)
--- a/ginac/inifcns_nstdsums.cpp
+++ b/ginac/inifcns_nstdsums.cpp
@@ -1043,6 +1043,25 @@ G_do_hoelder(std::vector<cln::cl_N> x, /* yes, it's passed by value */
        return result;
 }
 
+class less_object_for_cl_N
+{
+public:
+       bool operator() (const cln::cl_N & a, const cln::cl_N & b) const
+       {
+               // absolute value?
+               if (abs(a) != abs(b))
+                       return (abs(a) < abs(b)) ? true : false;
+
+               // complex phase?
+               if (phase(a) != phase(b))
+                       return (phase(a) < phase(b)) ? true : false;
+
+               // equal, therefore "less" is not true
+               return false;
+       }
+};
+
+
 // convergence transformation, used for numerical evaluation of G function.
 // the parameter x, s and y must only contain numerics
 static cln::cl_N
@@ -1050,17 +1069,17 @@ G_do_trafo(const std::vector<cln::cl_N>& x, const std::vector<int>& s,
            const cln::cl_N& y, bool flag_trailing_zeros_only)
 {
        // sort (|x|<->position) to determine indices
-       typedef std::multimap<cln::cl_R, std::size_t> sortmap_t;
+       typedef std::multimap<cln::cl_N, std::size_t, less_object_for_cl_N> sortmap_t;
        sortmap_t sortmap;
        std::size_t size = 0;
        for (std::size_t i = 0; i < x.size(); ++i) {
                if (!zerop(x[i])) {
-                       sortmap.insert(std::make_pair(abs(x[i]), i));
+                       sortmap.insert(std::make_pair(x[i], i));
                        ++size;
                }
        }
        // include upper limit (scale)
-       sortmap.insert(std::make_pair(abs(y), x.size()));
+       sortmap.insert(std::make_pair(y, x.size()));
 
        // generate missing dummy-symbols
        int i = 1;