From 3675da388f7ca12cdc6d93f4a8b94ef4bdd96b19 Mon Sep 17 00:00:00 2001
From: Alex Schmidt <alex@treefish.org>
Date: Wed, 26 Jun 2013 14:11:32 +0200
Subject: [PATCH] ...

---
 obstat.hpp | 35 ++++++++++++++++++++++++++++++++++-
 1 file changed, 34 insertions(+), 1 deletion(-)

diff --git a/obstat.hpp b/obstat.hpp
index 2996043..6818792 100644
--- a/obstat.hpp
+++ b/obstat.hpp
@@ -17,7 +17,9 @@ public:
 
   int computeEasy(const int& ival=0) { return mean(&measurements, ival); }
 
-  int computeJack(restype (*func)(vector< vector <meastype> > *vals, void *para)=NULL, void *para=NULL);
+  int computeJack(restype (*func)(vector< vector <meastype> > *vals, void *para), void *para=NULL);
+  int computeJack(void (*preMeasFunc)(vector< vector <meastype> > *allVals, vector <meastype> *preCalculated, void *para), 
+		  restype (*measFunc)(vector <meastype> *preCalculated, vector <meastype> *excludedmeas, void *para), void *para=NULL);
 
   restype getMean(int compid) { return computations[compid].val; }
   restype getErr(int compid) { return computations[compid].err; }
@@ -117,4 +119,35 @@ int obstat<meastype,restype>::computeJack(restype (*func)(vector< vector<meastyp
   return computations.size()-1;
 }
 
+template <typename meastype, typename restype>
+int obstat<meastype,restype>::computeJack(void (*preMeasFunc)(vector< vector <meastype> > *allVals, vector<meastype> *preCalculated, void *para), 
+					  restype (*measFunc)(vector<meastype> *preCalculated, vector<meastype> *excludedmeas, void *para), 
+					  void *para) {
+  int nmeas=measurements.size();
+  restype manymeans[nmeas];
+  result jackres;
+  vector<meastype> preCalculated;
+
+  jackres.val = 0;
+  jackres.err = 0;
+
+  preMeasFunc(&measurements, &preCalculated, para);
+
+  int imeas=0;
+  for(typename vector< vector<meastype> >::iterator removedIt = measurements.begin(); removedIt != measurements.end(); ++removedIt, imeas++) {
+    manymeans[imeas] = measFunc(&preCalculated, &(*removedIt), para);
+    jackres.val += manymeans[imeas]; 
+  }
+  jackres.val /= nmeas;
+  
+  for(int imean=0; imean<nmeas; imean++)
+    jackres.err += pow( manymeans[imean] - jackres.val, 2 );
+  jackres.err *= (double)(nmeas-1)/nmeas;
+  jackres.err = sqrt(jackres.err);
+
+  computations.push_back(jackres);
+
+  return computations.size()-1;
+}
+
 #endif
-- 
2.39.5