refactormain

1 files changed, 22 insertions, 23 deletions

diff --git a/ts.c b/ts.c
index a85bf0a..ad6081c 100644
--- a/ts.c
+++ b/ts.c
@@ -14,7 +14,7 @@
  * Uses a running accumulator to avoid floating-point drift,
  * and forces the last element to exactly x1.
  */
-fp_t* ts_fill(fp_t *x, uint_t n, fp_t x0, fp_t x1 )
+fp_t * ts_fill( fp_t *x, uint_t n, fp_t x0, fp_t x1 )
 {
 	if( !x || n == 0 ) return NULL;
 	if( n == 1 )
@@ -52,7 +52,6 @@ fp_t* ts_fill(fp_t *x, uint_t n, fp_t x0, fp_t x1 )
  */
 void ts_print( fp_t const *x, uint_t n, FILE *stream )
 {
-	assert( x != NULL );
 	if( !stream ) stream = stdout;
 	fputc( '[', stream );
 	if( !x || n == 0 )
@@ -61,13 +60,12 @@ void ts_print( fp_t const *x, uint_t n, FILE *stream )
 		return;
 	}
 
-	uint_t i;
-	for( i = 0; i < n - 1; i++ )
+	for( uint_t i = 0; i < n - 1; i++ )
 	{
 		fprintf( stream, "%.17g,", x[i] );
 	}
 
-	fprintf( stream, "%.17g]", x[i] );
+	fprintf( stream, "%.17g]", x[n - 1] );
 }
 
 
@@ -141,16 +139,15 @@ fp_t _ts_sum_to_covariance( fp_t sum_xy, fp_t sum_x, fp_t sum_y, uint_t k )
 {
 	if( k == 0 )
 		return 0.0;
-	fp_t result = ( sum_xy - ( sum_x * sum_y / k ) ) / k;
-	if( result < 0.0 )
-		result = 0.0;
-	return result;
+	return ( sum_xy - ( sum_x * sum_y / k ) ) / k;
 }
 
 
 fp_t _ts_sum_to_variance( fp_t sum_x2, fp_t sum_x, uint_t k )
 {
-	return _ts_sum_to_covariance( sum_x2, sum_x, sum_x, k );
+	/* Variance must be non-negative. Clamp to zero to handle
+	 * catastrophic cancellation in the naive summation formula. */
+	return fmax( 0.0, _ts_sum_to_covariance( sum_x2, sum_x, sum_x, k ) );
 }
 
 
@@ -165,7 +162,7 @@ uint_t _ts_circle_increment( uint_t i, uint_t k )
 /**
  * Allocate and zero-initialize a ts_constant_t structure.
  */
-ts_constant_t* ts_constant( void )
+ts_constant_t * ts_constant( void )
 {
 	return calloc( 1, sizeof( ts_constant_t ) );
 }
@@ -202,17 +199,23 @@ bool ts_test(
 	ts_constant_t const *c,
 	ts_parameter_t const *m )
 {
-	if ( !ts_function || !ts_offset )
+	if( !ts_function || !ts_offset )
 		return false;
 
 	if( !m ) return false;
 	if( !e ) return false;
 	if( !x ) return false;
 	if( m->n_constants > 0 && !c ) return false;
-	if( m->n_outputs > 1024 ) return false;
+	if( m->n_outputs > TS_MAX_LENGTH ) return false;
+
+	if( n > SIZE_MAX / sizeof( fp_t ) )
+		return false;
+
+	if( m->n_outputs > SIZE_MAX / sizeof( fp_t * ) )
+		return false;
 
 	uint_t offset = ts_offset( c );
-	if ( offset >= n )
+	if( offset >= n )
 	{
 		return false;
 	}
@@ -225,11 +228,6 @@ bool ts_test(
 		return false;
 	}
 
-	if ( n > UINT64_MAX / sizeof( fp_t ) )
-	{
-		return false;
-	}
-
 	for( uint_t i = 0; i < m->n_outputs; i++ )
 	{
 		fp_t *tmp = calloc( n, sizeof( fp_t ) );
@@ -260,7 +258,7 @@ bool ts_test(
 			fp_t ev = e[i][j];
 			if( isnan( yv ) || isnan( ev ) )
 			{
-				fprintf( stderr, "FAIL NaN y[%d][%d]=%f exp=%f\n", (int)i, (int)j, yv, ev );
+				fprintf( stderr, "FAIL NaN y[%" TS_PRIuint "][%" TS_PRIuint "]=%" TS_PRIfp " exp=%" TS_PRIfp "\n", i, j, yv, ev );
 				result = false;
 				continue;
 			}
@@ -269,7 +267,7 @@ bool ts_test(
 				continue;
 			if( isinf( yv ) || isinf( ev ) )
 			{
-				fprintf( stderr, "FAIL inf y[%d][%d]=%f exp=%f\n", (int)i, (int)j, yv, ev );
+				fprintf( stderr, "FAIL inf y[%" TS_PRIuint "][%" TS_PRIuint "]=%" TS_PRIfp " exp=%" TS_PRIfp "\n", i, j, yv, ev );
 				result = false;
 				continue;
 			}
@@ -277,9 +275,10 @@ bool ts_test(
 			fp_t abs_e = fabs( ev );
 			fp_t tol = TS_TEST_EPSILON + abs_e * TS_TEST_EPSILON;
 			if( diff > tol )
-				fprintf( stderr, "FAIL y[%d][%d]=%.6f exp=%.6f diff=%.8f tol=%.8f\n", (int)i, (int)j, yv, ev, diff, tol );
-			if( diff > tol )
+			{
+				fprintf( stderr, "FAIL y[%" TS_PRIuint "][%" TS_PRIuint "]=%.6f exp=%.6f diff=%.8f tol=%.8f\n", i, j, yv, ev, diff, tol );
 				result = false;
+			}
 		}
 
 		free( y[i] );