[Rivet-svn] rivet: Improve linspace (and hence also logspace) precision erro...

[Rivet Mercurial]

details:   https://rivet.hepforge.org/hg/rivet/rev/5d88921f711d
branches:  release-2-4-x
changeset: 4965:5d88921f711d
user:      Andy Buckley <andy at insectnation.org>
date:      Sun Dec 20 10:26:26 2015 +0000
description:
Improve linspace (and hence also logspace) precision errors by using multiplication rather than repeated addition to build edge list (thanks to Holger Schulz for the suggestion).

diffs (truncated from 51 to 50 lines):

--- a/ChangeLog	Thu Dec 17 21:40:37 2015 +0000
+++ b/ChangeLog	Sun Dec 20 10:26:26 2015 +0000
@@ -1,3 +1,9 @@
+2015-12-20  Andy Buckley  <andy.buckley at cern.ch>
+
+	* Improve linspace (and hence also logspace) precision errors by
+	using multiplication rather than repeated addition to build edge
+	list (thanks to Holger Schulz for the suggestion).
+
 2015-12-15  Andy Buckley  <andy.buckley at cern.ch>
 
 	* Add cmphistos and make-plots machinery for handling 'suffix'
--- a/include/Rivet/Math/MathUtils.hh	Thu Dec 17 21:40:37 2015 +0000
+++ b/include/Rivet/Math/MathUtils.hh	Sun Dec 20 10:26:26 2015 +0000
@@ -273,13 +273,11 @@
     assert(nbins > 0);
     vector<double> rtn;
     const double interval = (end-start)/static_cast<double>(nbins);
-    double edge = start;
     for (size_t i = 0; i < nbins; ++i) {
-      rtn.push_back(edge);
-      edge += interval;
+      rtn.push_back(start + i*interval);
     }
     assert(rtn.size() == nbins);
-    if (include_end) rtn.push_back(end); // exact end, not result of n * interval
+    if (include_end) rtn.push_back(end); //< exact end, not result of n * interval
     return rtn;
   }
 
@@ -295,15 +293,15 @@
     assert(nbins > 0);
     const double logstart = std::log(start);
     const double logend = std::log(end);
-    const vector<double> logvals = linspace(nbins, logstart, logend);
+    const vector<double> logvals = linspace(nbins, logstart, logend, false);
     assert(logvals.size() == nbins+1);
-    vector<double> rtn; rtn.reserve(logvals.size());
-    rtn.push_back(start);
-    for (size_t i = 1; i < logvals.size()-1; ++i) {
+    vector<double> rtn; rtn.reserve(nbins+1);
+    rtn.push_back(start); //< exact start, not exp(log(start))
+    for (size_t i = 1; i < logvals.size(); ++i) {
       rtn.push_back(std::exp(logvals[i]));
     }
     assert(rtn.size() == nbins);
-    if (include_end) rtn.push_back(end);
+    if (include_end) rtn.push_back(end); //< exact end, not exp(n * loginterval)
     return rtn;
   }