[PATCH] Run benchtest until given accuracy is reached

[OndÅej BÃlka <neleai at seznam dot cz>]

Hi, 

Benchtest now runs 10s which is often much more than necessary.

I used standard techique of calculating 
 http://en.wikipedia.org/wiki/Standard_error  
which with probability 95% bounds error caused by random 
factors only to 1%. I could tigthen these bounds if you wish.

This does not remove errors caused by biases which are hard to remove
by running benchmark longer. They must be carefully eliminated which is
future work.

OK to commit?

Ondra

---
 benchtests/Makefile         |    8 +--
 benchtests/bench-skeleton.c |  130 +++++++++++++++++++++++--------------------
 2 files changed, 72 insertions(+), 66 deletions(-)

diff --git a/benchtests/Makefile b/benchtests/Makefile
index 19e1be6..9e5d133 100644
--- a/benchtests/Makefile
+++ b/benchtests/Makefile
@@ -79,12 +79,8 @@ include ../Rules
 
 binaries-bench := $(addprefix $(objpfx)bench-,$(bench))
 
-# The default duration: 10 seconds.
-ifndef BENCH_DURATION
-BENCH_DURATION := 10
-endif
-
-CPPFLAGS-nonlib = -DDURATION=$(BENCH_DURATION)
+# We need sqrt to compute standard error.
+LDFLAGS += -lm
 
 # This makes sure CPPFLAGS-nonlib and CFLAGS-nonlib are passed
 # for all these modules.
diff --git a/benchtests/bench-skeleton.c b/benchtests/bench-skeleton.c
index 7359184..ee3ccf8 100644
--- a/benchtests/bench-skeleton.c
+++ b/benchtests/bench-skeleton.c
@@ -18,83 +18,93 @@
 
 #include <string.h>
 #include <stdint.h>
+#include <stdlib.h>
 #include <stdio.h>
 #include <time.h>
 #include <inttypes.h>
+#include <math.h>
 
-#define TIMESPEC_AFTER(a, b) \
-  (((a).tv_sec == (b).tv_sec) ?						      \
-     ((a).tv_nsec > (b).tv_nsec) :					      \
-	((a).tv_sec > (b).tv_sec))
 int
 main (int argc, char **argv)
 {
-  unsigned long i, k;
-  struct timespec start, end, runtime;
+  unsigned long k;
+  struct timespec start, end;
+	unsigned int seed = 42;
 
-  memset (&runtime, 0, sizeof (runtime));
   memset (&start, 0, sizeof (start));
   memset (&end, 0, sizeof (end));
 
   clock_getres (CLOCK_PROCESS_CPUTIME_ID, &start);
 
-  /* Measure 1000 times the resolution of the clock.  So for a 1ns resolution
-     clock, we measure 1000 iterations of the function call at a time.
-     Measurements close to the minimum clock resolution won't make much sense,
-     but it's better than having nothing at all.  */
-  unsigned long iters = 1000 * start.tv_nsec;
+	/* We do measuement repeately to make clock overhead small.  */
+  unsigned long iters = 10 * start.tv_nsec;
+  unsigned long total_iters = 1000;
+  double d_total_i;
+	double mean, variance, standard_error;
+  uint64_t total, max, min;
+  int variant;
 
-  for (int v = 0; v < NUM_VARIANTS; v++)
+  for (variant = 0; variant < NUM_VARIANTS; variant++)
     {
-      /* Run for approximately DURATION seconds.  */
-      clock_gettime (CLOCK_MONOTONIC_RAW, &runtime);
-      runtime.tv_sec += DURATION;
-
-      double d_total_i = 0;
-      uint64_t total = 0, max = 0, min = 0x7fffffffffffffff;
-      while (1)
-	{
-	  for (i = 0; i < NUM_SAMPLES (v); i++)
-	    {
-	      clock_gettime (CLOCK_PROCESS_CPUTIME_ID, &start);
-	      for (k = 0; k < iters; k++)
-		BENCH_FUNC (v, i);
-	      clock_gettime (CLOCK_PROCESS_CPUTIME_ID, &end);
-
-	      uint64_t cur = (end.tv_nsec - start.tv_nsec
-			      + ((end.tv_sec - start.tv_sec)
-				 * (uint64_t) 1000000000));
-
-	      if (cur > max)
-		max = cur;
-
-	      if (cur < min)
-		min = cur;
-
-	      total += cur;
-
-	      d_total_i += iters;
-	    }
-	  struct timespec curtime;
-
-	  memset (&curtime, 0, sizeof (curtime));
-	  clock_gettime (CLOCK_MONOTONIC_RAW, &curtime);
-	  if (TIMESPEC_AFTER (curtime, runtime))
-	    goto done;
-	}
-
-      double d_total_s;
-      double d_iters;
-
-    done:
-      d_total_s = total * 1e-9;
-      d_iters = iters;
+      /* Double number of iterations until we reach 1% precision with
+         95% confidence.  */
+      do
+        {
+          total_iters *= 2;
+          mean = 0;
+          double delta = 0;
+          double M2 = 0;
+          d_total_i = 0;
+          total = 0;
+          max = 0;
+          min = 0x7fffffffffffffff;
+          long n;
+          for (n = 1; n <= total_iters; n++)
+            {
+                  clock_gettime (CLOCK_PROCESS_CPUTIME_ID, &start);
+                  int sample = rand_r (&seed) % (NUM_SAMPLES (variant));
+                  for (k = 0; k < iters; k++)
+                    BENCH_FUNC (variant, sample);
+                  clock_gettime (CLOCK_PROCESS_CPUTIME_ID, &end);
+
+                  uint64_t cur = (end.tv_nsec - start.tv_nsec
+                                  + ((end.tv_sec - start.tv_sec)
+                                     * (uint64_t) 1000000000));
+
+                  /* Numericaly stable computation of variance.  See:
+                     Donald E. Knuth (1998) The Art of Computer Programming,
+                     volume 2: Seminumerical Algorithms, 3rd edn., p. 232.
+                     Boston: Addison-Wesley.  */
+                  delta = cur - mean;
+                  mean = mean + delta / n;
+                  M2 = M2 + delta * (cur - mean);
+
+                  if (cur > max)
+                    max = cur;
+
+                  if (cur < min)
+                    min = cur;
+
+                  total += cur;
+
+                  d_total_i += iters;
+            }
+          /* See: http://en.wikipedia.org/wiki/Variance  */
+          variance = M2 / (n - 1);
+          /* See: http://en.wikipedia.org/wiki/Standard_error  */
+          standard_error = sqrt (variance) / sqrt(n);
+          /* As we double iterations probability of false positive
+             decreases exponentialy with factor sqrt(2).  */
+        }
+      while (1.96 * standard_error * (sqrt(2) - 1) > mean / 100.0);
+
+      double d_total_s = total * 1e-9;
+      double d_iters = iters;
 
       printf ("%s: ITERS:%g: TOTAL:%gs, MAX:%gns, MIN:%gns, %g iter/s\n",
-	      VARIANT (v),
-	      d_total_i, d_total_s, max / d_iters, min / d_iters,
-	      d_total_i / d_total_s);
+              VARIANT (variant),
+              d_total_i, d_total_s, max / d_iters, min / d_iters,
+              d_total_i / d_total_s);
     }
-
   return 0;
 }
-- 
1.7.4.4