1 files changed, 29 insertions, 24 deletions

diff --git a/ts.h b/ts.h
index fdb522e..43761ea 100644
--- a/ts.h
+++ b/ts.h
@@ -232,8 +232,7 @@ static inline int_t _ts_read_int_t( const uint8_t *p )
 
 /*
  * ts_fill — Linearly fill array x with values from x0 to x1.
- * RETURNS: the same pointer x (for convenience). Caller retains ownership.
- * PRECONDITION: x must not be NULL; n > 0.
+ * RETURNS: x on success, NULL if x is NULL or n == 0. Caller retains ownership.
  */
 fp_t * ts_fill( fp_t *x, uint_t n, fp_t x0, fp_t x1 );
 
@@ -278,44 +277,50 @@ fp_t _ts_sum_to_variance( fp_t, fp_t, uint_t );
  * PRECONDITION: pointer w must not be NULL for all ts_welford_* functions.
  */
 typedef struct {
-    fp_t mean;
-    fp_t m2;
-    uint_t count;
+	fp_t mean;
+	fp_t m2;
+	uint_t count;
 } ts_welford_t;
 
 static inline void ts_welford_init( ts_welford_t *w )
 {
-    if (!w) return;
-    w->mean = 0; w->m2 = 0; w->count = 0;
+	if ( !w ) return;
+	w->mean = 0; w->m2 = 0; w->count = 0;
 }
 
 static inline void ts_welford_add( ts_welford_t *w, fp_t x )
 {
-    if (!w) return;
-    w->count++;
-    fp_t delta = x - w->mean;
-    w->mean += delta / w->count;
-    fp_t delta2 = x - w->mean;
-    w->m2 += delta * delta2;
+	if ( !w ) return;
+	w->count++;
+	fp_t delta = x - w->mean;
+	w->mean += delta / w->count;
+	fp_t delta2 = x - w->mean;
+	w->m2 += delta * delta2;
 }
 
 static inline void ts_welford_remove( ts_welford_t *w, fp_t x )
 {
-    if (!w) return;
-    if (w->count <= 1) { ts_welford_init(w); return; }
-    fp_t old_mean = w->mean;
-    w->count--;
-    w->mean = (old_mean * (w->count + 1) - x) / w->count;
-    fp_t delta = x - old_mean;
-    fp_t delta2 = x - w->mean;
-    w->m2 -= delta * delta2;
-    if (w->m2 < 0) w->m2 = 0;
+	if ( !w ) return;
+	if ( w->count <= 1 ) { ts_welford_init( w ); return; }
+	fp_t old_mean = w->mean;
+	w->count--;
+	w->mean = old_mean + ( old_mean - x ) / w->count;
+	fp_t delta = x - old_mean;
+	fp_t delta2 = x - w->mean;
+	w->m2 -= delta * delta2;
+	w->m2 = fmax( 0.0, w->m2 );
 }
 
 static inline fp_t ts_welford_variance( ts_welford_t *w )
 {
-    if (!w) return 0;
-    return w->count > 0 ? w->m2 / w->count : 0;
+	if ( !w ) return 0;
+	return w->count > 0 ? w->m2 / w->count : 0;
+}
+
+static inline fp_t ts_welford_sample_variance( ts_welford_t *w )
+{
+	if ( !w ) return 0;
+	return w->count > 1 ? w->m2 / ( w->count - 1 ) : 0;
 }