diff --git a/Contributors.md b/Contributors.md
new file mode 100644
index 0000000..7d36357
--- /dev/null
+++ b/Contributors.md
@@ -0,0 +1 @@
+bliksemstraal
\ No newline at end of file
diff --git a/stats.go b/stats.go
index ab22363..b9685b7 100644
--- a/stats.go
+++ b/stats.go
@@ -10,107 +10,107 @@ import (
 // fitness of genomes in the population.
 type Stats struct {
 	max, min float64
-	mean     float64
-	sumsq    float64 // sum of squares of deviation from the mean
+	sum      float64
+	sumsq    float64
 	count    float64
 }
 
 // Put inserts a new value into the data.
-func (s Stats) Put(x float64) Stats {
-	if s.count == 0 {
-		s.max = math.Inf(-1)
-		s.min = math.Inf(+1)
+func (s *Stats) Put(x float64) Stats {
+	if x > s.max || s.count == 0 {
+		s.max = x
+	}
+	if x < s.min || s.count == 0 {
+		s.min = x
 	}
 
-	delta := x - s.mean
-	newcount := s.count + 1
-
-	// max & min
-	s.max = math.Max(s.max, x)
-	s.min = math.Min(s.min, x)
-
-	// mean
-	s.mean += delta / newcount
-
-	// sum of squares
-	s.sumsq += delta * delta * (s.count / newcount)
-
-	// count
-	s.count = newcount
+	s.sum += x
+	s.sumsq += x * x
+	s.count++
 
-	return s
+	return *s
 }
 
 // Merge merges the data of two Stats objects.
-func (s Stats) Merge(t Stats) Stats {
-	if s.count == 0 {
-		s.max = math.Inf(-1)
-		s.min = math.Inf(+1)
-	}
-
-	delta := t.mean - s.mean
-	newcount := t.count + s.count
-
-	// max & min
-	s.max = math.Max(s.max, t.max)
-	s.min = math.Min(s.min, t.min)
+func (s *Stats) Merge(t Stats) Stats {
 
-	// mean
-	s.mean += delta * (t.count / newcount)
+	// Do not merge if t is empty.
+	if t.count == 0 {
+		return *s
+	}
 
-	// sum of squares
+	s.count += t.count
+	s.sum += t.sum
 	s.sumsq += t.sumsq
-	s.sumsq += delta * delta * (t.count * s.count / newcount)
 
-	// count
-	s.count = newcount
+	// The OR clause is used in case s is empty (max and min will be 0)
+	// and t.max < 0 or t.min > 0.
+	if t.max > s.max || s.count == 0 {
+		s.max = t.max
+	}
+	if t.min > s.min || s.count == 0 {
+		s.min = t.min
+	}
 
-	return s
+	return *s
 }
 
 // Max returns the maximum data point.
-func (s Stats) Max() float64 {
+func (s *Stats) Max() float64 {
 	return s.max
 }
 
 // Min returns the minimum data point.
-func (s Stats) Min() float64 {
+func (s *Stats) Min() float64 {
 	return s.min
 }
 
 // Range returns the difference in the maximum and minimum data points.
-func (s Stats) Range() float64 {
+func (s *Stats) Range() float64 {
 	return s.max - s.min
 }
 
 // Mean returns the average of the data.
-func (s Stats) Mean() float64 {
-	return s.mean
+func (s *Stats) Mean() float64 {
+	if s.count == 0 {
+		return 0
+	}
+	return s.sum / s.count
 }
 
-// Var returns the population variance of the data.
-func (s Stats) Var() float64 {
-	return s.sumsq / s.count
+// Var returns the population variance of the data calculated in a single-
+// pass method.
+func (s *Stats) Var() float64 {
+	if s.count < 2 {
+		return 0 // no variation
+	}
+	return (s.count*s.sumsq - s.sum*s.sum) / (s.count * s.count)
 }
 
 // SD returns the population standard deviation of the data.
-func (s Stats) SD() float64 {
-	return math.Sqrt(s.sumsq / s.count)
+func (s *Stats) SD() float64 {
+	if s.count < 2 {
+		return 0 // no variation
+	}
+	return math.Sqrt(s.Var())
 }
 
 // RSD returns the population relative standard deviation of the data, also
 // known as the coefficient of variation.
-func (s Stats) RSD() float64 {
+func (s *Stats) RSD() float64 {
+	if s.count == 0 {
+		return 0
+	}
 	return s.SD() / s.Mean()
 }
 
 // Count returns the size of the data.
-func (s Stats) Count() int {
+func (s *Stats) Count() int {
 	return int(s.count)
 }
 
 // String returns a string listing a summary of the statistics.
-func (s Stats) String() string {
+func (s *Stats) String() string {
 	return fmt.Sprintf("Max: %f | Min: %f | SD: %f",
 		s.Max(),
 		s.Min(),