G_do_hoelder: fix bug related to comparing a float with a rational number in CLN.

[ginac.git] / ginac / inifcns_nstdsums.cpp

diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp
index 474ed3e233b79adcd94235fd09ca703d902dc7b6..c6f5012afbab8da2305854cec7094d340b8540a9 100644 (file)
--- a/ginac/inifcns_nstdsums.cpp
+++ b/ginac/inifcns_nstdsums.cpp
@@ -47,7 +47,7 @@
  */
 
 /*
- *  GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2021 Johannes Gutenberg University Mainz, Germany
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -1013,21 +1013,25 @@ G_do_hoelder(std::vector<cln::cl_N> x, /* yes, it's passed by value */
        for (std::size_t i = 0; i < size; ++i)
                x[i] = x[i]/y;
 
+        // 24.03.2021: this block can be outside the loop over r 
+       cln::cl_RA p(2);
+       bool adjustp;
+       do {
+               adjustp = false;
+               for (std::size_t i = 0; i < size; ++i) {
+                        // 24.03.2021: replaced (x[i] == cln::cl_RA(1)/p) by (cln::zerop(x[i] - cln::cl_RA(1)/p)
+                        //             in the case where we compare a float with a rational, CLN behaves differently in the two versions
+                       if (cln::zerop(x[i] - cln::cl_RA(1)/p) ) {
+                               p = p/2 + cln::cl_RA(3)/2;
+                               adjustp = true;
+                               continue;
+                       }
+               }
+       } while (adjustp);
+       cln::cl_RA q = p/(p-1);
+
        for (std::size_t r = 0; r <= size; ++r) {
                cln::cl_N buffer(1 & r ? -1 : 1);
-               cln::cl_RA p(2);
-               bool adjustp;
-               do {
-                       adjustp = false;
-                       for (std::size_t i = 0; i < size; ++i) {
-                               if (x[i] == cln::cl_RA(1)/p) {
-                                       p = p/2 + cln::cl_RA(3)/2;
-                                       adjustp = true;
-                                       continue;
-                               }
-                       }
-               } while (adjustp);
-               cln::cl_RA q = p/(p-1);
                std::vector<cln::cl_N> qlstx;
                std::vector<int> qlsts;
                for (std::size_t j = r; j >= 1; --j) {