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.
+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
// convergence transformation, used for numerical evaluation of G function.
// the parameter x, s and y must only contain numerics
static cln::cl_N
const cln::cl_N& y, bool flag_trailing_zeros_only)
{
// sort (|x|<->position) to determine indices
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_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)
++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;
// generate missing dummy-symbols
int i = 1;