diff --git a/ChangeLog b/ChangeLog
index 6470cb51..f6c3e50e 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -26,6 +26,8 @@
 - im_remainder() has a float result for float types
 - im_measure() allows sel == NULL, meaning all patches
 - added "deprecated" package
+- faster, simpler, better im_max()
+- im_max() returns true modulus, not square of modulus, for complex images
 
 25/3/09 started 7.18.0
 - revised version numbers
diff --git a/TODO b/TODO
index 1c6e0aa2..f5a27593 100644
--- a/TODO
+++ b/TODO
@@ -1,4 +1,17 @@
-- im_stats() needs a rewrite ... scrap the stupid 42 thing
+- new im_max() needs testing
+
+  im_avg() needs to work for complex: return modulus
+
+  use im_max() model to rewrite im_min(), im_stats()
+
+- im_minpos(), im_maxpos() need partialing
+
+- im_maxpos_vec needs gtkdoc
+
+
+
+
+
 
 - 1-bit PNG read is broken?
 
@@ -10,29 +23,6 @@
 
 
 
-- im__cast_and_call() no longer does
-
-	out->Bbits = im_bits_of_fmt( out->BandFmt );
-
-  make sure we do this for each use of im__cast_and_call()
-
-- im__cast_and_call() no longer casts input to output type, it casts inputs 
-  up to the smallest common type
-  
-  so we need to change all the buffer processors from 
-  
-  	T op T -> T
-
-  to 
-
-  	T1 op T1 -> T2
-
-  done add/sub/mul, doing div
-  
-  others?
-
-
-
 - rename vipsCC in SWIG as pyvips?
 
 
diff --git a/libvips/arithmetic/im_avg.c b/libvips/arithmetic/im_avg.c
index 04629ab6..1b735645 100644
--- a/libvips/arithmetic/im_avg.c
+++ b/libvips/arithmetic/im_avg.c
@@ -100,6 +100,19 @@ stop_fn( void *seq, void *a, void *b )
 	return( 0 );
 }
 
+/* Sum pels in this section.
+ */
+#define LOOP( TYPE ) { \
+	TYPE *p; \
+	\
+	for( y = to; y < bo; y++ ) { \
+		p = (TYPE *) IM_REGION_ADDR( reg, le, y ); \
+		\
+		for( x = 0; x < sz; x++ ) \
+			sum += *p++; \
+	} \
+}
+
 /* Loop over region, accumulating a sum in *tmp.
  */
 static int
@@ -115,19 +128,6 @@ scan_fn( REGION *reg, void *seq, void *a, void *b )
 	double sum = 0.0;
 	int x, y;
 
-/* Sum pels in this section.
- */
-#define LOOP( TYPE ) { \
-	TYPE *p; \
-	\
-	for( y = to; y < bo; y++ ) { \
-		p = (TYPE *) IM_REGION_ADDR( reg, le, y ); \
-		\
-		for( x = 0; x < sz; x++ ) \
-			sum += *p++; \
-	} \
-}
-
 	/* Now generate code for all types. 
 	 */
 	switch( im->BandFmt ) {
diff --git a/libvips/arithmetic/im_linreg.c b/libvips/arithmetic/im_linreg.c
index bf6aba26..95bcc5f8 100644
--- a/libvips/arithmetic/im_linreg.c
+++ b/libvips/arithmetic/im_linreg.c
@@ -1,12 +1,4 @@
-/* @(#) Function to find perform pixelwise linear regression on an array of 
- * @(#) single band images.
- * @(#)
- * @(#) int im_linreg(
- * @(#)   IMAGE **ins,
- * @(#)   IMAGE *out,
- * @(#)   double *xs
- * @(#) );
- * @(#) 
+/* im_linreg.c
  *
  * Copyright: 2006, The Nottingham Trent University
  *
@@ -16,6 +8,8 @@
  *
  * 12/5/09
  *	- make x_anal() static, fix some signed/unsigned warnings
+ * 3/9/09
+ * 	- gtkdoc comment
  */
 
 /*
@@ -148,6 +142,16 @@ SKIP_ALL_DECL( double );
 
 /** EXPORTED FUNCTION DEFINITION **/
 
+/** im_linreg:
+ * @ins: NULL-terminated array of input images
+ * @out: results of analysis
+ * @xs:	X position of each image (pixel value is Y)
+ *
+ * Function to find perform pixelwise linear regression on an array of 
+ * single band images. The output is a seven-band douuble image
+ *
+ * TODO: figure out how this works and fix up these docs!
+ */
 int im_linreg( IMAGE **ins, IMAGE *out, double *xs ){
 #define FUNCTION_NAME "im_linreg"
   int n;
diff --git a/libvips/arithmetic/im_max.c b/libvips/arithmetic/im_max.c
index 2713432d..81e8b8bd 100644
--- a/libvips/arithmetic/im_max.c
+++ b/libvips/arithmetic/im_max.c
@@ -18,6 +18,9 @@
  *	- random wrong result for >1 thread :-( (thanks Joe)
  * 15/10/07
  * 	- oh, heh, seq->inf was not being set correctly, not that it mattered
+ * 4/9/09
+ * 	- rewrite with im__value(), much simpler and fixes a race condition
+ * 	- gtkdoc comment
  */
 
 /*
@@ -58,7 +61,6 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
-#include <assert.h>
 
 #include <vips/vips.h>
 
@@ -66,39 +68,15 @@
 #include <dmalloc.h>
 #endif /*WITH_DMALLOC*/
 
-/* Per-call state.
- */
-typedef struct _MaxInfo {
-	/* Parameters.
-	 */
-	IMAGE *in;
-	double *out;
-
-	/* Global max so far.
-	 */
-	double value;
-	int valid;		/* zero means value is unset */
-} MaxInfo;
-
-/* Per thread state.
- */
-typedef struct _Seq {
-	MaxInfo *inf;
-
-	double value;
-	int valid;		/* zero means value is unset */
-} Seq;
-
 /* New sequence value.
  */
 static void *
 max_start( IMAGE *in, void *a, void *b )
 {
-	MaxInfo *inf = (MaxInfo *) a;
-	Seq *seq = IM_NEW( NULL, Seq );
+	double *value = (double *) a;
+	double *seq = IM_NEW( NULL, double );
 
-	seq->inf = inf;
-	seq->valid = 0;
+	*seq = *value;
 
 	return( (void *) seq );
 }
@@ -108,47 +86,19 @@ max_start( IMAGE *in, void *a, void *b )
 static int
 max_stop( void *vseq, void *a, void *b )
 {
-	Seq *seq = (Seq *) vseq;
-	MaxInfo *inf = (MaxInfo *) a;
+	double *seq = (double *) vseq;
+	double *value = (double *) a;
 
-	if( seq->valid ) {
-		if( !inf->valid )
-			/* Just copy.
-			 */
-			inf->value = seq->value;
-		else 
-			/* Merge.
-			 */
-			inf->value = IM_MAX( inf->value, seq->value );
-
-		inf->valid = 1;
-	}
+	/* Merge.
+	 */
+	*value = IM_MAX( *value, *seq );
 
 	im_free( seq );
 
 	return( 0 );
 }
 
-/* Loop over region, adding to seq.
- */
-static int
-max_scan( REGION *reg, void *vseq, void *a, void *b )
-{
-	Seq *seq = (Seq *) vseq;
-	Rect *r = &reg->valid;
-	IMAGE *im = reg->im;
-	int le = r->left;
-	int to = r->top;
-	int bo = IM_RECT_BOTTOM(r);
-	int nel = IM_REGION_N_ELEMENTS( reg );
-
-	int x, y;
-
-	double m;
-
-#define loop(TYPE) { \
-	m = *((TYPE *) IM_REGION_ADDR( reg, le, to )); \
-	\
+#define LOOP( TYPE ) { \
 	for( y = to; y < bo; y++ ) { \
 		TYPE *p = (TYPE *) IM_REGION_ADDR( reg, le, y ); \
 		\
@@ -161,22 +111,17 @@ max_scan( REGION *reg, void *vseq, void *a, void *b )
 	} \
 } 
 
-#define complex_loop(TYPE) { \
-	TYPE *p = (TYPE *) IM_REGION_ADDR( reg, le, to ); \
-	double real = p[0]; \
-	double imag = p[1]; \
-	\
-	m = real * real + imag * imag; \
-	\
+#define CLOOP( TYPE ) { \
 	for( y = to; y < bo; y++ ) { \
 		TYPE *p = (TYPE *) IM_REGION_ADDR( reg, le, y ); \
 		\
-		for( x = 0; x < nel * 2; x += 2 ) { \
-			double mod; \
+		for( x = 0; x < nel; x++ ) { \
+			double mod, re, im; \
 			\
-			real = p[x]; \
-			imag = p[x + 1]; \
-			mod = real * real + imag * imag; \
+			re = p[0]; \
+			im = p[1]; \
+			p += 2; \
+			mod = re * re + im * im; \
 			\
 			if( mod > m ) \
 				m = mod; \
@@ -184,47 +129,76 @@ max_scan( REGION *reg, void *vseq, void *a, void *b )
 	} \
 } 
 
+/* Loop over region, adding to seq.
+ */
+static int
+max_scan( REGION *reg, void *vseq, void *a, void *b )
+{
+	double *seq = (double *) vseq;
+	Rect *r = &reg->valid;
+	IMAGE *im = reg->im;
+	int le = r->left;
+	int to = r->top;
+	int bo = IM_RECT_BOTTOM(r);
+	int nel = IM_REGION_N_ELEMENTS( reg );
+
+	int x, y;
+	double m;
+
+	m = *seq;
+
 	switch( im->BandFmt ) {
-	case IM_BANDFMT_UCHAR:		loop( unsigned char ); break; 
-	case IM_BANDFMT_CHAR:		loop( signed char ); break; 
-	case IM_BANDFMT_USHORT:		loop( unsigned short ); break; 
-	case IM_BANDFMT_SHORT:		loop( signed short ); break; 
-	case IM_BANDFMT_UINT:		loop( unsigned int ); break;
-	case IM_BANDFMT_INT:		loop( signed int ); break; 
-	case IM_BANDFMT_FLOAT:		loop( float ); break; 
-	case IM_BANDFMT_DOUBLE:		loop( double ); break; 
-	case IM_BANDFMT_COMPLEX:	complex_loop( float ); break; 
-	case IM_BANDFMT_DPCOMPLEX:	complex_loop( double ); break; 
+	case IM_BANDFMT_UCHAR:		LOOP( unsigned char ); break; 
+	case IM_BANDFMT_CHAR:		LOOP( signed char ); break; 
+	case IM_BANDFMT_USHORT:		LOOP( unsigned short ); break; 
+	case IM_BANDFMT_SHORT:		LOOP( signed short ); break; 
+	case IM_BANDFMT_UINT:		LOOP( unsigned int ); break;
+	case IM_BANDFMT_INT:		LOOP( signed int ); break; 
+	case IM_BANDFMT_FLOAT:		LOOP( float ); break; 
+	case IM_BANDFMT_DOUBLE:		LOOP( double ); break; 
+	case IM_BANDFMT_COMPLEX:	CLOOP( float ); break; 
+	case IM_BANDFMT_DPCOMPLEX:	CLOOP( double ); break; 
 
 	default:  
-		assert( 0 );
+		g_assert( 0 );
 	}
 
-	if( seq->valid ) {
-		seq->value = IM_MAX( seq->value, m ); 
-	}
-	else {
-		seq->value = m;
-		seq->valid = 1;
-	}
+	*seq = m; 
 
 #ifdef DEBUG
         printf( "im_max: left = %d, top = %d, width = %d, height = %d\n",
 		r->left, r->top, r->width, r->height );
-	printf( "   (max = %g)\n", seq->value );
+	printf( "   (max = %g)\n", *seq );
 #endif /*DEBUG*/
 
 	return( 0 );
 }
 
+/* Get the value of pixel (0, 0). Use this to init the min/max value for
+ * threads. Shared with im_min(), im_stats() etc. This will return mod for
+ * complex.
+ */
+int
+im__value( IMAGE *im, double *value )
+{
+	IMAGE *t;
+
+	if( !(t = im_open_local( im, "im__value", "p" )) ||
+		im_extract_area( im, t, 0, 0, 1, 1 ) ||
+		im_avg( t, value ) )
+		return( -1 );
+
+	return( 0 );
+}
+
 /** 
  * im_max:
  * @in: input #IMAGE
  * @out: output double
  *
  * Finds the the maximum value of image #in and returns it at the
- * location pointed by out.  If input is complex, the max square amplitude 
- * (re*re+im*im) is returned. im_max() finds the maximum of all bands: if you
+ * location pointed by out.  If input is complex, the max modulus
+ * is returned. im_max() finds the maximum of all bands: if you
  * want to find the maximum of each band separately, use im_stats().
  *
  * See also: im_maxpos(), im_min(), im_stats().
@@ -233,21 +207,29 @@ max_scan( REGION *reg, void *vseq, void *a, void *b )
  */
 int
 im_max( IMAGE *in, double *out )
-{	
-	MaxInfo inf;
-
-	inf.in = in;
-	inf.out = out;
-	inf.valid = 0;
+{
+	double value;
 
 	if( im_pincheck( in ) ||
 		im_check_uncoded( "im_max", in ) )
 		return( -1 );
 
-	if( im_iterate( in, max_start, max_scan, max_stop, &inf, NULL ) ) 
+	if( im__value( in, &value ) )
+		return( -1 );
+	/* We use square mod for scanning, for speed.
+	 */
+	if( im_iscomplex( in ) )
+		value *= value;
+
+	if( im_iterate( in, max_start, max_scan, max_stop, &value, NULL ) ) 
 		return( -1 );
 
-	*out = inf.value;
+	/* Back to modulus.
+	 */
+	if( im_iscomplex( in ) )
+		value = sqrt( value );
+
+	*out = value;
 
 	return( 0 );
 }
diff --git a/libvips/arithmetic/im_maxpos.c b/libvips/arithmetic/im_maxpos.c
index 3a25214f..dec838e8 100644
--- a/libvips/arithmetic/im_maxpos.c
+++ b/libvips/arithmetic/im_maxpos.c
@@ -1,15 +1,4 @@
-/* @(#) Function to find the maximum of an image. Works for any 
- * @(#) image type. Returns a double and the location of max.
- * @(#)
- * @(#) Function im_maxpos() assumes that input
- * @(#) is either memory mapped or in a buffer.
- * @(#)
- * @(#) int im_maxpos(in, xpos, ypos, max)
- * @(#) IMAGE *in;
- * @(#) int *xpos, *ypos;
- * @(#) double *max;
- * @(#)
- * @(#) Returns 0 on success and -1 on error
+/* im_maxpos.c
  *
  * Copyright: 1990, J. Cupitt
  *
@@ -21,6 +10,8 @@
  *	- im_incheck() call added
  * 20/6/95 JC
  *	- now returns double for value, like im_max()
+ * 4/9/09
+ * 	- gtkdoc comment
  */
 
 /*
@@ -55,7 +46,6 @@
 #include <vips/intl.h>
 
 #include <stdio.h>
-#include <assert.h>
 
 #include <vips/vips.h>
 
@@ -77,8 +67,22 @@
 	case IM_BANDFMT_COMPLEX:	loopcmplx(float); break; \
 	case IM_BANDFMT_DPCOMPLEX:	loopcmplx(double); break; 
 
-/* Find the position of the maximum of an image. Take any format, returns a 
- * float and its position.
+/**
+ * im_maxpos:
+ * @in: image to search
+ * @xpos: returned x position of maximum
+ * @ypos: returned y position of maximum
+ * @out: returned pixel value at that position
+ *
+ * Function to find the maximum of an image. Works for any 
+ * image type. Returns a double and the location of max.
+ *
+ * This is not a PIO operation! It may use a lot of memory and take a while.
+ * Needs a rewrite.
+ *
+ * See also: im_minpos(), im_min(), im_stats(), im_maxpos_avg().
+ *
+ * Returns: 0 on success, -1 on error
  */
 int
 im_maxpos( IMAGE *in, int *xpos, int *ypos, double *out )
@@ -135,7 +139,7 @@ im_maxpos( IMAGE *in, int *xpos, int *ypos, double *out )
 	switch( in->BandFmt ) {
 		im_for_all_types();
 		default: { 
-			assert( 0 );
+			g_assert( 0 );
 			return( -1 );
 		}
 	}
diff --git a/libvips/arithmetic/im_min.c b/libvips/arithmetic/im_min.c
index 55a76c70..75cb095f 100644
--- a/libvips/arithmetic/im_min.c
+++ b/libvips/arithmetic/im_min.c
@@ -1,11 +1,4 @@
-/* @(#) Function to find the minimum of an image. Works for any 
- * @(#) image type. Returns a double.
- * @(#)
- * @(#) int im_min(in, min)
- * @(#) IMAGE *in;
- * @(#) double *min;
- * @(#)
- * @(#) Returns 0 on success and -1 on error
+/* im_min.c
  *
  * Copyright: 1990, J. Cupitt
  *
@@ -221,6 +214,20 @@ scan_fn( REGION *reg, void *vseq, void *a, void *b )
 	return( 0 );
 }
 
+/** 
+ * im_min:
+ * @in: input #IMAGE
+ * @out: output double
+ *
+ * Finds the the minimum value of image @in and returns it at the
+ * location pointed by @out.  If input is complex, the min square amplitude 
+ * (re*re+im*im) is returned. im_min() finds the maximum of all bands: if you
+ * want to find the maximum of each band separately, use im_stats().
+ *
+ * See also: im_minpos(), im_max(), im_stats().
+ *
+ * Returns: 0 on success, -1 on error
+ */
 int
 im_min( IMAGE *in, double *out )
 {	
@@ -230,12 +237,9 @@ im_min( IMAGE *in, double *out )
 	inf.out = out;
 	inf.valid = 0;
 
-	if( im_pincheck( in ) )
+	if( im_pincheck( in ) ||
+		im_check_uncoded( "im_min", in ) )
 		return( -1 );
-	if( in->Coding != IM_CODING_NONE ) {
-		im_error( "im_min", "%s", _( "not uncoded" ) );
-		return( -1 );
-	}
 
 	if( im_iterate( in, start_fn, scan_fn, stop_fn, &inf, NULL ) ) 
 		return( -1 );
diff --git a/libvips/arithmetic/im_minpos.c b/libvips/arithmetic/im_minpos.c
index 6861d69a..72ccc3e5 100644
--- a/libvips/arithmetic/im_minpos.c
+++ b/libvips/arithmetic/im_minpos.c
@@ -1,15 +1,4 @@
-/* @(#) Function to find the minimim of an image. Works for any 
- * @(#) image type. Returns a double and the location of min
- * @(#)
- * @(#) Function im_minpos() assumes that input
- * @(#) is either memory mapped or in a buffer.
- * @(#)
- * @(#) int im_minpos(in, xpos, ypos, out)
- * @(#) IMAGE *in;
- * @(#) int *xpos, *ypos;
- * @(#) double *out;
- * @(#)
- * @(#) Returns 0 on success and -1 on error
+/* im_minpos.c
  *
  * Copyright: 1990, J. Cupitt
  *
@@ -22,6 +11,8 @@
  *	- im_incheck() added
  * 20/6/95 JC
  *	- now returns double for value, like im_max()
+ * 4/9/09
+ * 	- gtkdoc comment
  */
 
 /*
@@ -56,7 +47,6 @@
 #include <vips/intl.h>
 
 #include <stdio.h>
-#include <assert.h>
 
 #include <vips/vips.h>
 
@@ -64,7 +54,23 @@
 #include <dmalloc.h>
 #endif /*WITH_DMALLOC*/
 
-/* Find the minimum of an image. Take any format, returns a double.  */
+/**
+ * im_minpos:
+ * @in: image to search
+ * @xpos: returned x position of minimum
+ * @ypos: returned y position of minimum
+ * @out: returned pixel value at that position
+ *
+ * Function to find the minimum of an image. Works for any 
+ * image type. Returns a double and the location of min.
+ *
+ * This is not a PIO operation! It may use a lot of memory and take a while.
+ * Needs a rewrite.
+ *
+ * See also: im_maxpos(), im_max(), im_stats(), im_maxpos_avg().
+ *
+ * Returns: 0 on success, -1 on error
+ */
 int
 im_minpos( IMAGE *in, int *xpos, int *ypos, double *out )
 {
@@ -131,7 +137,7 @@ im_minpos( IMAGE *in, int *xpos, int *ypos, double *out )
 	case IM_BANDFMT_DPCOMPLEX:	loopcmplx(double); break; 
 
 	default: 
-		assert( 0 );
+		g_assert( 0 );
 	}
 
 	/* Take out bands on x.
diff --git a/libvips/include/vips/internal.h b/libvips/include/vips/internal.h
index 1321de1d..44041e34 100644
--- a/libvips/include/vips/internal.h
+++ b/libvips/include/vips/internal.h
@@ -135,6 +135,7 @@ int im__cast_and_call( IMAGE *in1, IMAGE *in2, IMAGE *out,
 	im_wrapmany_fn fn, void *a );
 VipsBandFmt im__format_common( IMAGE *in1, IMAGE *in2 );
 int im__math( const char *name, IMAGE *in, IMAGE *out, im_wrapone_fn gen );
+int im__value( IMAGE *im, double *value );
 
 int im__test_kill( IMAGE *im );
 void *im__mmap( int fd, int writeable, size_t length, gint64 offset );