initial import into cvs

source is JT personal cvs ~/src/libmisc++/


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@31 f88db872-0e4f-0410-b76b-b9085cfa78c5

6 files changed, 971 insertions, 0 deletions

diff --git a/src/jtutil/cpm_map.cc b/src/jtutil/cpm_map.cc
new file mode 100644
index 0000000..db6ac8e
--- /dev/null
+++ b/src/jtutil/cpm_map.cc
@@ -0,0 +1,99 @@
+// cpm_map.cc -- "integer +/-" mapping  i --> j = const +/- i
+// $Id$
+//
+// cpm_map::{plus,minus}
+// cpm_map::cpm_map	// mirror map, specified by fixed point
+// cpm_map::cpm_map	// generic map specified by sample point & sign
+// cpm_map::cpm_map	// generic map specified by *double* sample point & sign
+//
+
+#include <assert.h>
+#include <stdio.h>
+#include <jt/stdc.h>
+#include <jt/util++.hh>
+#include <jt/cpm_map.hh>
+
+//*****************************************************************************
+//*****************************************************************************
+// cpm_map -- ipm mapping from integers <--> floating point values
+//*****************************************************************************
+//*****************************************************************************
+
+//
+// we declared and initilized these constants in the .hh file, but
+// we still have to *define* them here
+//
+const bool cpm_map::plus_map;
+const bool cpm_map::minus_map;
+
+//*****************************************************************************
+
+//
+// This function constructs a  cpm_map  object with a "-" sign and a
+// specified fixed point (must be integer or half-integer) and domain.
+//
+cpm_map::cpm_map(int min_i_in, int max_i_in,
+		 double fixed_point)
+	: min_i_(min_i_in), max_i_(max_i_in),
+	  sign_(minus_map)
+{
+double d_offset = 2.0 * fixed_point;
+if (! fuzzy<double>::is_integer(d_offset))
+   then error_exit(ERROR_EXIT,
+"***** cpm_map::cpm_map (mirror):\n"
+"        fixed_point=%g isn't (fuzzily) integral or half-integral!\n"
+,
+		   fixed_point);				/*NOTREACHED*/
+
+offset_ = round<double>::to_integer(d_offset);
+}
+
+//*****************************************************************************
+
+//
+// This function constructs a generic  cpm_map  object, with the mapping
+// specified by a sample point  sample_i --> sample_j  and by sign (one
+// of  {plus,minus}_map ).  The sample point need not be in domain/range.
+//
+cpm_map::cpm_map(int min_i_in, int max_i_in,
+		 int sample_i, int sample_j,
+		 bool sign_in)
+		 : min_i_(min_i_in), max_i_(max_i_in),
+		   offset_(sign_in ? sample_j - sample_i
+				   : sample_j + sample_i),
+		   sign_(sign_in)
+{
+// should be guaranteed by above  offset_  setup
+assert( map(sample_i) == sample_j );
+}
+
+//*****************************************************************************
+
+//
+// This function constructs a generic  cpm_map  object, with the mapping
+// specified by a *double* sample point  sample_i --> sample_j  (which
+// must specify an  integer --> integer  mapping, i.e.  4.2 --> 4.2  is
+// ok for a + map, and 4.5 --> 4.5 is ok for a minus map, but  4.2 --> 4.7
+// is never ok) and by sign (one of  {plus,minus}_map ).  The sample point
+// need not be in domain/range.
+//
+cpm_map::cpm_map(int min_i_in, int max_i_in,
+		 double sample_i, double sample_j,
+		 bool sign_in)
+		 : min_i_(min_i_in), max_i_(max_i_in),
+		   sign_(sign_in)
+{
+double d_offset = sign_in ? sample_j - sample_i
+			  : sample_j + sample_i;
+if (! fuzzy<double>::is_integer(d_offset))
+   then error_exit(ERROR_EXIT,
+"***** cpm_map::cpm_map (generic via fp sample point):\n"
+"        d_offset=%g isn't fuzzily integral!\n"
+"        ==> sample_i=%g --> sample_j=%g\n"
+"            doesn't fuzzily specify an  integer --> integer  mapping!\n"
+,
+		   d_offset,
+		   sample_i, sample_j);				/*NOTREACHED*/
+
+offset_ = round<double>::to_integer(d_offset);
+}
diff --git a/src/jtutil/cpm_map.hh b/src/jtutil/cpm_map.hh
new file mode 100644
index 0000000..ee09fb8
--- /dev/null
+++ b/src/jtutil/cpm_map.hh
@@ -0,0 +1,121 @@
+// cpm_map.hh -- "integer +/-" mapping  i --> j = const +/- i
+// $Id$
+//
+// cpm_map - "integer +/-" mapping
+//
+
+//
+// prerequisites
+//    <assert.h>
+//
+
+//
+// This class represents a unit-slope linear integer mapping
+//	i --> j = const +/- i
+//
+
+#ifndef NDEBUG
+//
+// Full bounds checking is done on the mapping in both directions.
+//
+#endif
+
+//*****************************************************************************
+
+class	cpm_map
+	{
+public:
+	// bounds info -- domain
+	int min_i() const { return min_i_; }
+	int max_i() const { return max_i_; }
+	int N_points() const { return HOW_MANY_IN_RANGE(min_i_,max_i_); }
+	bool in_domain(int i) const { return i >= min_i() && i <= max_i(); }
+
+	// is the mapping + or - ?
+	static const bool plus_map  = true;
+	static const bool minus_map = false;
+	bool sign_flag() const { return sign_; }
+	int sign() const { return sign_ ? +1 : -1; }
+
+	// the mapping itself
+	int map_unchecked(int i) const
+		{
+		return sign_flag() ? offset_ + i
+				   : offset_ - i;
+		}
+	int inv_map_unchecked(int j) const
+		{
+		return sign_flag() ? j - offset_
+				   : offset_ - j;
+		}
+	int map(int i) const
+		{
+		assert(in_domain(i));
+		return map_unchecked(i);
+		}
+	int inv_map(int j) const
+		{
+		int i = inv_map_unchecked(j);
+		assert(in_domain(i));
+		return i;
+		}
+
+	// bounds info -- range
+	int min_j() const
+		{
+		return sign_flag() ? map_unchecked(min_i())
+				   : map_unchecked(max_i());
+		}
+	int max_j() const
+		{
+		return sign_flag() ? map_unchecked(max_i())
+				   : map_unchecked(min_i());
+		}
+	bool in_range(int j) const { return in_domain(inv_map_unchecked(j)); }
+
+	//
+	// constructors
+	//
+
+	// "mirror" map: i --> const - i
+	// ... map specified by fixed point (must be integer or half-integer)
+	cpm_map(int min_i_in, int max_i_in,
+		double fixed_point);
+
+	// "shift" map: i --> const + i
+	// ... map specified by shift amount
+	// ... default is identity map
+	cpm_map(int min_i_in, int max_i_in,
+		int shift_amount = 0)
+		: min_i_(min_i_in), max_i_(max_i_in),
+		  offset_(shift_amount), sign_(plus_map)
+		{ }
+
+	// generic map: i --> const +/- i
+	// ... map specified by sample point sample_i --> sample_j
+	//     and by sign (one of  {plus,minus}_map )
+	// ... sample point is *not* checked for being in domain/range
+	cpm_map(int min_i_in, int max_i_in,
+		int sample_i, int sample_j,
+		bool sign_in);
+
+	// generic map: i --> const +/- i
+	// ... map specified by *double* sample point sample_i --> sample_j
+	//     (must specify an integer --> integer mapping)
+	//     and by sign (one of  {plus,minus}_map )
+	// ... sample point is *not* checked for being in domain/range
+	cpm_map(int min_i_in, int max_i_in,
+		double sample_i, double sample_j,
+		bool sign_in);
+
+	// no need for explicit destructor, compiler-generated no-op is ok
+	// ditto for copy constructor and assignment operator
+
+private:
+	// bounds (inclusive)
+	int min_i_, max_i_;
+
+	// these define the actual mapping
+	int offset_;
+	bool sign_;
+	};
diff --git a/src/jtutil/linear_map.cc b/src/jtutil/linear_map.cc
new file mode 100644
index 0000000..518c3c3
--- /dev/null
+++ b/src/jtutil/linear_map.cc
@@ -0,0 +1,249 @@
+// linear_map.cc -- linear mapping from integers <--> floating point values
+// $Id$
+//
+// linear_map::linear_map - basic constructor
+// linear_map::linear_map - constructor for subrange of existing linear_map
+// linear_map::fp_int_of_fp - convert fp --> int coord, return as fp
+// linear_map::int_of_fp - convert fp --> int, check for being fuzzily integral
+// linear_map::delta_int_of_delta_fp - ... same for spacings (delta_* values)
+//
+// *** instantiation of linear_map template for <double>
+//
+
+#include <assert.h>
+#include <stdio.h>
+#include <jt/stdc.h>
+#include <jt/util++.hh>
+#include <jt/linear_map.hh>
+
+//*****************************************************************************
+//*****************************************************************************
+// linear_map -- linear mapping from integers <--> floating point values
+//*****************************************************************************
+//*****************************************************************************
+
+//
+// This function constructs a  linear_map<fp>  object.
+//
+template <class fp>
+linear_map<fp>::linear_map(int min_int_in, int max_int_in,
+			   fp min_fp_in, fp delta_fp_in, fp max_fp_in)
+	: delta_(delta_fp_in), inverse_delta_(1.0 / delta_fp_in),
+	  min_int_(min_int_in), max_int_(max_int_in)
+{
+constructor_common(min_fp_in, max_fp_in);
+}
+
+//*****************************************************************************
+
+//
+// This function constructs a  linear_map<fp>  object with a subrange
+// of an existing one.
+//
+template <class fp>
+linear_map<fp>::linear_map(const linear_map<fp> &lm_in,
+			   int min_int_in, int max_int_in)	// subrange
+	: delta_(lm_in.delta_fp()), inverse_delta_(lm_in.inverse_delta_fp()),
+	  min_int_(min_int_in), max_int_(max_int_in)
+{
+if (! (in_range(min_int_in) && in_range(max_int_in)) )
+   then error_exit(ERROR_EXIT,
+"***** linear_map<fp>::linear_map:\n"
+"        min_int_in=%d and/or max_int_in=%d\n"
+"        aren't in integer range [%d,%d] of existing linear_map!\n"
+,
+		   min_int_, max_int_,
+		   lm_in.min_int(), lm_in.max_int());		/*NOTREACHED*/
+
+constructor_common(lm_in.fp_of_int_unchecked(min_int_in),
+		   lm_in.fp_of_int_unchecked(max_int_in));
+}
+
+//*****************************************************************************
+
+//
+// This function does the common argument validation and setup for
+// all the constructors of class  linear_map<fp>:: .
+//
+template <class fp>
+void linear_map<fp>::constructor_common(fp min_fp_in, fp max_fp_in)
+	// assumes min_int_, max_int_, delta_, inverse_delta_ already initialized
+	// ==> ok to use min_int(), max_int(), delta_fp(), inverse_delta_fp()
+	// ... other class members *not* yet initialized
+{
+offset_ = 0.0;		// temp value
+offset_ = min_fp_in - fp_of_int_unchecked(min_int());
+
+// this should be guaranteed by the above calculation
+assert(fuzzy<fp>::EQ(fp_of_int_unchecked(min_int()), min_fp_in));
+
+// this is a test of the consistency of the input arguments
+if (fuzzy<fp>::NE(fp_of_int_unchecked(max_int()), max_fp_in))
+   then error_exit(ERROR_EXIT,
+"***** linear_map<fp>::linear_map:\n"
+"        int range [%d,%d]\n"
+"        and fp range [%g(%g)%g]\n"
+"        are (fuzzily) inconsistent!\n"
+,
+		   min_int(), max_int(),
+		   double(min_fp_in), double(delta_fp()), double(max_fp_in));
+								/*NOTREACHED*/
+}
+
+//*****************************************************************************
+
+//
+// This function converts  fp  --> int  coordinate, returning the result
+// as an fp (which need not be fuzzily integral).
+//
+template <class fp>
+fp linear_map<fp>::fp_int_of_fp(fp x)
+	const
+{
+if (! in_range(x))
+   then error_exit(ERROR_EXIT,
+"***** linear_map<fp>::fp_int_of_fp:\n"
+"        fp value x=%g is (fuzzily) outside the grid!\n"
+"        {min(delta)max}_fp = %g(%g)%g\n"
+,
+		   double(x),
+		   double(min_fp()), double(delta_fp()), double(max_fp()));
+								/*NOTREACHED*/
+
+return inverse_delta_ * (x - offset_);
+}
+
+//*****************************************************************************
+
+//
+// This function converts  fp  --> int  and checks that the result is
+// fuzzily integral.  (The  nia  argument specifies what to do if the
+// result *isn't* fuzzily integral.)
+//
+template <class fp>
+int linear_map<fp>::int_of_fp(fp x, noninteger_action nia = error)
+	const
+{
+const fp fp_int = fp_int_of_fp(x);
+
+if (fuzzy<fp>::is_integer(fp_int))
+   then {
+	// x is (fuzzily) a grid point ==> return that
+	return ::round<fp>::to_integer(fp_int);		// *** EARLY RETURN ***
+	}
+
+// get to here ==> x isn't (fuzzily) a grid point
+static const char *const noninteger_msg =
+   "%s linear_map<fp>::int_of_fp:\n"
+   "        x=%g isn't (fuzzily) a grid point!\n"
+   "        {min(delta)max}_fp() = %g(%g)%g\n"
+   ;
+switch	(nia)
+	{
+case error:
+	error_exit(ERROR_EXIT,
+		   noninteger_msg,
+		   "*****",
+		   double(x),
+		   double(min_fp()), double(delta_fp()), double(max_fp()));
+								/*NOTREACHED*/
+
+case warning:
+	printf(noninteger_msg,
+	       "---",
+	       double(x),
+	       double(min_fp()), double(delta_fp()), double(max_fp()));
+	// fall through
+
+case round:
+	return ::round<fp>::to_integer(fp_int);		// *** EARLY RETURN ***
+
+case floor:
+	return ::round<fp>::floor(fp_int);		// *** EARLY RETURN ***
+
+case ceiling:
+	return ::round<fp>::ceiling(fp_int);		// *** EARLY RETURN ***
+
+default:
+	error_exit(PANIC_EXIT,
+"***** linear_map<fp>::int_of_fp: illegal nia=(int)%d\n"
+"                                 (this should never happen!)\n"
+,
+		   int(nia));					/*NOTREACHED*/
+	}
+return 0;	// dummy return to quiet gcc
+		// (which doesn't grok that error_exit() never returns)
+}
+
+//*****************************************************************************
+
+//
+// This function converts "delta" spacings in the fp coordinate to
+// corresponding "delta" spacings in the int coordinate, and checks that
+// the result is fuzzily integral.  (The  nia  argument specifies what to
+// do if the result *isn't* fuzzily integral.)
+//
+template <class fp>
+int linear_map<fp>::delta_int_of_delta_fp(fp delta_x,
+					  noninteger_action nia = error)
+	const
+{
+const fp fp_delta_int = inverse_delta_ * delta_x;
+
+if (fuzzy<fp>::is_integer(fp_delta_int))
+   then {
+	// delta_x is (fuzzily) an integer number of grid spacings
+	// ==> return that
+	return ::round<fp>::to_integer(fp_delta_int);	// *** EARLY RETURN ***
+	}
+
+// get to here ==> delta_x isn't (fuzzily) an integer number of grid spacings
+static const char *const noninteger_msg =
+   "%s linear_map<fp>::delta_int_of_delta_fp:\n"
+   "        delta_x=%g isn't (fuzzily) an integer number of grid spacings!\n"
+   "        {min(delta)max}_fp() = %g(%g)%g\n"
+   ;
+switch	(nia)
+	{
+case error:
+	error_exit(ERROR_EXIT,
+		   noninteger_msg,
+		   "*****",
+		   double(delta_x),
+		   double(min_fp()), double(delta_fp()), double(max_fp()));
+								/*NOTREACHED*/
+
+case warning:
+	printf(noninteger_msg,
+	       "---",
+	       double(delta_x),
+	       double(min_fp()), double(delta_fp()), double(max_fp()));
+	// fall through
+
+case round:
+	return ::round<fp>::to_integer(fp_delta_int);	// *** EARLY RETURN ***
+
+case floor:
+	return ::round<fp>::floor(fp_delta_int);	// *** EARLY RETURN ***
+
+case ceiling:
+	return ::round<fp>::ceiling(fp_delta_int);	// *** EARLY RETURN ***
+
+default:
+	error_exit(PANIC_EXIT,
+"***** linear_map<fp>::delta_int_of_delta_fp: illegal nia=(int)%d\n"
+"                                 (this should never happen!)\n"
+,
+		   int(nia));					/*NOTREACHED*/
+	}
+return 0;	// dummy return to quiet gcc
+		// (which doesn't grok that error_exit() never returns)
+}
+
+//*****************************************************************************
+//*****************************************************************************
+// *** instantiation of linear_map template for <double>
+//*****************************************************************************
+//*****************************************************************************
+
+template class linear_map<double>;
diff --git a/src/jtutil/linear_map.hh b/src/jtutil/linear_map.hh
new file mode 100644
index 0000000..f7ae2a8
--- /dev/null
+++ b/src/jtutil/linear_map.hh
@@ -0,0 +1,132 @@
+// linear_map.hh -- linear mapping from integers <--> floating point values
+// $Id$
+//
+// linear_map - linear mapping from integers <--> floating point values
+//
+
+//
+// prerequisites
+//    <assert.h>
+//    <jt/util++.hh>	// for fuzzy::
+//
+
+//
+// The template defined in this file represents a linear mapping between
+// a contiguous range of integers, and an arithmetic progression of
+// floating point values, parameterized by the floating point data type.
+//
+
+#ifndef NDEBUG
+//
+// Full bounds checking is done on the mapping in both directions
+//
+#endif
+
+//*****************************************************************************
+
+template <class fp>
+class	linear_map
+	{
+public:
+	// integer bounds info
+	int min_int() const { return min_int_; }
+	int max_int() const { return max_int_; }
+	int N_points() const { return HOW_MANY_IN_RANGE(min_int_,max_int_); }
+	bool in_range(int i) const
+		{ return (i >= min_int()) && (i <= max_int()); }
+	int clamp(int i) const
+		{
+		if	(i < min_int())
+		   then return min_int();
+		else if (i > max_int())
+		   then return max_int();
+		else	return i;
+		}
+
+	// convert int --> fp
+	fp fp_of_int_unchecked(int i) const
+		{ return offset_ + delta_*i; }
+	fp fp_of_int(int i) const
+		{
+		assert(in_range(i));
+		return fp_of_int_unchecked(i);
+		}
+
+	// converg delta_int --> delta_fp
+	fp delta_fp_of_delta_int(int delta_i) const
+		{ return delta_ * delta_i; }
+
+	// fp bounds info
+	fp min_fp() const { return fp_of_int_unchecked(min_int_); }
+	fp delta_fp() const { return delta_; }
+	fp inverse_delta_fp() const { return inverse_delta_; }
+	fp max_fp() const { return fp_of_int_unchecked(max_int_); }
+	bool in_range(fp x) const
+		{
+		return fuzzy<fp>::GE(x,min_fp()) && fuzzy<fp>::LE(x,max_fp());
+		}
+	fp clamp(fp x) const
+		{
+		if	(x < min_fp())
+		   then return min_fp();
+		else if (x > max_fp())
+		   then return max_fp();
+		else	return x;
+		}
+
+	// convert linear map indices <--> C-style 0-origin indices
+	int zero_origin_int(int i) const { return i - min_int(); }
+	int map_int(int zero_origin_i) { return zero_origin_i + min_int(); }
+
+
+
+	// convert fp --> int coordinate, but return result as fp
+	// (which need not be fuzzily integral)
+	fp fp_int_of_fp(fp x) const;
+
+	// convert fp --> int, check being fuzzily integral
+	enum	noninteger_action	// what to do if `int'
+					// isn't fuzzily integral?
+		{
+		error,		// error_exit(...)
+		warning,	// print warning msg, then round to nearest
+		round,		// (silently) round to nearest
+		floor,		// (silently) round to -infinity
+		ceiling		// (silently) round to +infinity
+		};
+	int int_of_fp(fp x, noninteger_action nia = error) const;
+
+	// convert delta_fp --> delta_int, check being fuzzily integral
+	int delta_int_of_delta_fp(fp delta_x, noninteger_action nia = error)
+		const;
+
+	// constructors
+	linear_map(int min_int_in, int max_int_in,
+		   fp min_fp_in, fp delta_fp_in, fp max_fp_in);
+	// ... construct with subrange of existing linear_map
+	linear_map(const linear_map<fp> &lm_in,
+		   int min_int_in, int max_int_in);
+
+	// no need for explicit destructor, compiler-generated no-op is ok
+
+	// no need for copy constructor or assignment operator,
+	// compiler-generated defaults are ok
+
+private:
+	// common code (argument validation & setup) for all constructors
+	// assumes min_int_, max_int_, delta_ already initialized,
+	//         other class members *not* initialized
+	void constructor_common(fp min_fp_in, fp max_fp_in);
+
+	// these define the actual mapping
+	// via the  fp_of_int()  function (above)
+	fp offset_, delta_;
+
+	// cache of 1.0/delta_
+	// ==> avoids fp divide in inverse_delta_fp()
+	// ==> also makes fp --> int conversions slightly faster
+	fp inverse_delta_;
+
+	// bounds (inclusive)
+	const int min_int_, max_int_;
+	};
diff --git a/src/jtutil/test_cpm_map.cc b/src/jtutil/test_cpm_map.cc
new file mode 100644
index 0000000..e4e0711
--- /dev/null
+++ b/src/jtutil/test_cpm_map.cc
@@ -0,0 +1,213 @@
+// test_cpm_map.cc -- test driver for cpm_map class
+// $Id$
+
+#include <assert.h>
+#include <stdio.h>
+#include <jt/stdc.h>
+#include <jt/util++.hh>
+#include <jt/cpm_map.hh>
+
+//*****************************************************************************
+
+// prototypes
+
+void test_shift_map(int min_i, int sample_i, int max_i,
+		    int min_j, int sample_j, int max_j);
+void test_mirror_map(int min_i, int sample_i, int max_i,
+		     int min_j, int sample_i, int max_j,
+		     double fixed_point);
+
+void verify_shift_map(const cpm_map sm,
+		      const int min_i, const int sample_i, const int max_i,
+		      const int min_j, const int sample_j, const int max_j);
+void verify_mirror_map(const cpm_map mm,
+		       const int min_i, const int sample_i, const int max_i,
+		       const int min_j, const int sample_j, const int max_j,
+		       const double fixed_point);
+
+//*****************************************************************************
+//*****************************************************************************
+//*****************************************************************************
+
+int main()
+{
+test_shift_map(4 , 10, 14,
+	       42, 48, 52);
+
+test_mirror_map(4 , 10, 14,
+		16, 20, 26,
+		15.0);
+test_mirror_map(4 , 10, 14,
+		17, 21, 27,
+		15.5);
+
+return 0;
+}
+
+//*****************************************************************************
+
+//
+// this function tests a shift map
+//
+void test_shift_map(const int min_i, const int sample_i, const int max_i,
+		    const int min_j, const int sample_j, const int max_j)
+{
+printf("testing shift_map [%d,%d,%d] --> [%d,%d,%d]...\n",
+       min_i, sample_i, max_i,
+       min_j, sample_j, max_j);
+cpm_map sm(min_i, max_i,
+	   sample_j - sample_i);
+verify_shift_map(sm,
+		 min_i, sample_i, max_i,
+		 min_j, sample_j, max_j);
+
+{
+printf("testing generic_map (+) [%d,%d,%d] --> [%d,%d,%d]...\n",
+       min_i, sample_i, max_i,
+       min_j, sample_j, max_j);
+cpm_map gm(min_i, max_i,
+	   sample_i, sample_j,
+	   cpm_map::plus_map);
+verify_shift_map(gm,
+		 min_i, sample_i, max_i,
+		 min_j, sample_j, max_j);
+}
+
+{
+const double delta = 0.42;
+const double d_sample_i = double(sample_i) + delta;
+const double d_sample_j = double(sample_j) + delta;
+printf("testing generic_map (+) [%d,%g,%d] --> [%d,%g,%d]...\n",
+       min_i, d_sample_i, max_i,
+       min_j, d_sample_j, max_j);
+cpm_map gm(min_i, max_i,
+	   d_sample_i, d_sample_j,
+	   cpm_map::plus_map);
+verify_shift_map(gm,
+		 min_i, sample_i, max_i,
+		 min_j, sample_j, max_j);
+}
+}
+
+//*****************************************************************************
+
+//
+// this function tests a mirror map
+//
+void test_mirror_map(const int min_i, const int sample_i, const int max_i,
+		     const int min_j, const int sample_j, const int max_j,
+		     const double fixed_point)
+{
+{
+printf("testing mirror_map [%d,%d,%d] --> [%d,%d,%d] with fixed point %g\n",
+       min_i, sample_i, max_i,
+       min_j, sample_j, max_j,
+       fixed_point);
+cpm_map mm(min_i, max_i,
+	   fixed_point);
+verify_mirror_map(mm,
+		  min_i, sample_i, max_i,
+		  min_j, sample_j, max_j,
+		  fixed_point);
+}
+
+{
+printf("testing generic_map (-) [%d,%d,%d] --> [%d,%d,%d]...\n",
+       min_i, sample_i, max_i,
+       min_j, sample_j, max_j);
+cpm_map gm(min_i, max_i,
+	   sample_i, sample_j,
+	   cpm_map::minus_map);
+verify_mirror_map(gm,
+		  min_i, sample_i, max_i,
+		  min_j, sample_j, max_j,
+		  fixed_point);
+}
+
+{
+const double delta = 0.5;
+const double d_sample_i = double(sample_i) + delta;
+const double d_sample_j = double(sample_j) - delta;
+printf("testing generic_map (-) [%d,%g,%d] --> [%d,%g,%d]...\n",
+       min_i, d_sample_i, max_i,
+       min_j, d_sample_j, max_j);
+cpm_map gm(min_i, max_i,
+	   d_sample_i, d_sample_j,
+	   cpm_map::minus_map);
+verify_mirror_map(gm,
+		  min_i, sample_i, max_i,
+		  min_j, sample_j, max_j,
+		  fixed_point);
+}
+}
+
+//*****************************************************************************
+//*****************************************************************************
+//*****************************************************************************
+
+//
+// this function verifies that a shift map is correct
+//
+void verify_shift_map(const cpm_map sm,
+		      const int min_i, const int sample_i, const int max_i,
+		      const int min_j, const int sample_j, const int max_j)
+{
+assert(sm.min_i() == min_i);		assert(sm.min_j() == min_j);
+assert(sm.max_i() == max_i);		assert(sm.max_j() == max_j);
+assert(sm.N_points() == HOW_MANY_IN_RANGE(min_i, max_i));
+
+assert(! sm.in_domain(min_i-1));	assert(! sm.in_range(min_j-1));
+assert(  sm.in_domain(min_i  ));	assert(  sm.in_range(min_j  ));
+assert(  sm.in_domain(min_i+1));	assert(  sm.in_range(min_j+1));
+assert(  sm.in_domain(max_i-1));	assert(  sm.in_range(max_j-1));
+assert(  sm.in_domain(max_i  ));	assert(  sm.in_range(max_j  ));
+assert(! sm.in_domain(max_i+1));	assert(! sm.in_range(max_j+1));
+
+assert(sm.sign_flag() == cpm_map::plus_map);
+assert(sm.sign() == +1);
+
+assert(sm.map(min_i) == min_j);
+assert(sm.map(max_i) == max_j);
+assert(sm.map(sample_i) == sample_j);
+
+assert(sm.inv_map(min_j) == min_i);
+assert(sm.inv_map(max_j) == max_i);
+assert(sm.inv_map(sample_j) == sample_i);
+
+printf("==> all ok!\n");
+}
+
+//*****************************************************************************
+
+//
+// this function verifies that a mirror map is correct
+//
+void verify_mirror_map(const cpm_map mm,
+		       const int min_i, const int sample_i, const int max_i,
+		       const int min_j, const int sample_j, const int max_j,
+		       const double fixed_point)
+{
+assert(mm.min_i() == min_i);		assert(mm.min_j() == min_j);
+assert(mm.max_i() == max_i);		assert(mm.max_j() == max_j);
+assert(mm.N_points() == HOW_MANY_IN_RANGE(min_i, max_i));
+
+assert(! mm.in_domain(min_i-1));	assert(! mm.in_range(min_j-1));
+assert(  mm.in_domain(min_i  ));	assert(  mm.in_range(min_j  ));
+assert(  mm.in_domain(min_i+1));	assert(  mm.in_range(min_j+1));
+assert(  mm.in_domain(max_i-1));	assert(  mm.in_range(max_j-1));
+assert(  mm.in_domain(max_i  ));	assert(  mm.in_range(max_j  ));
+assert(! mm.in_domain(max_i+1));	assert(! mm.in_range(max_j+1));
+
+assert(mm.sign_flag() == cpm_map::minus_map);
+assert(mm.sign() == -1);
+
+assert(mm.map(min_i) == max_j);
+assert(mm.map(max_i) == min_j);
+assert(mm.map(sample_i) == sample_j);
+
+assert(mm.inv_map(min_j) == max_i);
+assert(mm.inv_map(max_j) == min_i);
+assert(mm.inv_map(sample_j) == sample_i);
+
+printf("==> all ok!\n");
+}
diff --git a/src/jtutil/test_linear_map.cc b/src/jtutil/test_linear_map.cc
new file mode 100644
index 0000000..1e1dc73
--- /dev/null
+++ b/src/jtutil/test_linear_map.cc
@@ -0,0 +1,157 @@
+// test_linear_map.cc -- test driver for linear_map class
+// $Id$
+
+#include <assert.h>
+#include <stdio.h>
+#include <jt/stdc.h>
+#include <jt/util++.hh>
+#include <jt/linear_map.hh>
+
+//*****************************************************************************
+
+//
+// Usage:
+//	test_linear_map
+// If an argument is supplied (it doesn't matter what it is), then
+// the array indices will be printed as they're tested.
+//
+int main(int argc, const char *const argv[])
+{
+const int min_int = 4;		const int max_int = 14;
+const int middle_int = (min_int + max_int) / 2;
+const int min_widen = 13;	const int max_widen = 9;
+const double min_fp = 9.0;
+const double delta_fp = 0.1;
+const double max_fp = 10.0;
+const double middle_fp = (min_fp + max_fp) / 2.0;
+
+printf("constructing linear_map<double>...\n");
+linear_map<double> lm_wide(min_int-min_widen, max_int+max_widen,
+			   min_fp-min_widen*delta_fp,
+			   delta_fp,
+			   max_fp+max_widen*delta_fp);
+linear_map<double> lm = *new linear_map<double>(lm_wide, min_int, max_int);
+
+printf("checking bounds and range info...\n");
+assert(lm.min_int() == min_int);
+assert(lm.max_int() == max_int);
+assert(lm.N_points() == HOW_MANY_IN_RANGE(min_int,max_int));
+
+assert(   lm.in_range(min_int)   );
+assert(   lm.in_range(max_int)   );
+assert( ! lm.in_range(min_int-1) );
+assert( ! lm.in_range(max_int+1) );
+
+assert( lm.clamp(min_int-100) == min_int    );
+assert( lm.clamp(min_int-1  ) == min_int    );
+assert( lm.clamp(min_int    ) == min_int    );
+assert( lm.clamp(min_int+1  ) == min_int+1  );
+assert( lm.clamp(middle_int ) == middle_int );
+assert( lm.clamp(max_int-1  ) == max_int-1  );
+assert( lm.clamp(max_int    ) == max_int    );
+assert( lm.clamp(max_int+1  ) == max_int    );
+assert( lm.clamp(max_int+100) == max_int    );
+
+assert( fuzzy<double>::EQ(lm.min_fp(), min_fp) );
+assert( fuzzy<double>::EQ(lm.delta_fp(), delta_fp) );
+assert( fuzzy<double>::EQ(lm.inverse_delta_fp(), 1.0/delta_fp) );
+assert( fuzzy<double>::EQ(lm.max_fp(), max_fp) );
+
+assert(   lm.in_range(min_fp)         );
+assert(   lm.in_range(max_fp)         );
+assert( ! lm.in_range(min_fp-1.0e-10) );
+assert( ! lm.in_range(max_fp+1.0e-10) );
+
+assert( fuzzy<double>::EQ(lm.clamp(min_fp - 100.0*delta_fp), min_fp) );
+assert( fuzzy<double>::EQ(lm.clamp(min_fp -   0.4*delta_fp), min_fp) );
+assert( fuzzy<double>::EQ(lm.clamp(min_fp), min_fp) );
+assert( fuzzy<double>::EQ(lm.clamp(min_fp + 0.4*delta_fp),
+				   min_fp + 0.4*delta_fp) );
+assert( fuzzy<double>::EQ(lm.clamp(middle_fp), middle_fp) );
+assert( fuzzy<double>::EQ(lm.clamp(max_fp - 0.4*delta_fp),
+				   max_fp - 0.4*delta_fp) );
+assert( fuzzy<double>::EQ(lm.clamp(max_fp), max_fp) );
+assert( fuzzy<double>::EQ(lm.clamp(max_fp +   0.4*delta_fp), max_fp) );
+assert( fuzzy<double>::EQ(lm.clamp(max_fp + 100.0*delta_fp), max_fp) );
+
+printf("checking delta roundings...\n");
+    for (int i = 1 ; i <= 3 ; ++i)
+    {
+    printf("   testing i = %d...\n", i);
+    assert( fuzzy<double>::EQ( lm.delta_fp_of_delta_int(i), i*delta_fp ) );
+    assert( lm.delta_int_of_delta_fp(i*delta_fp) == i );
+    double di = double(i);
+
+    // test silent roundings
+    assert( lm.delta_int_of_delta_fp((di-0.49)*delta_fp, lm.round) == i );
+    assert( lm.delta_int_of_delta_fp((di+0.49)*delta_fp, lm.round) == i );
+
+    // test floor/ceiling
+    assert( lm.delta_int_of_delta_fp((di-0.01)*delta_fp, lm.floor  ) == i-1 );
+    assert( lm.delta_int_of_delta_fp((di-0.01)*delta_fp, lm.ceiling) == i   );
+    assert( lm.delta_int_of_delta_fp((di+0.01)*delta_fp, lm.round  ) == i   );
+    assert( lm.delta_int_of_delta_fp((di+0.01)*delta_fp, lm.ceiling) == i+1 );
+    }
+
+printf("this should produce a pair of warning messages...\n");
+assert( lm.delta_int_of_delta_fp(2.99*delta_fp, lm.warning) == 3 );
+assert( lm.delta_int_of_delta_fp(3.01*delta_fp, lm.warning) == 3 );
+
+printf("checking mappings...\n");
+int i;		// must declare these outside loop
+double x;	// because we have *two* loop variables!
+	for (i = min_int, x = min_fp ;
+	     i <= max_int ;
+	     ++i, x += delta_fp)
+	{
+	printf("   testing i = %d of [%d,%d]...\n", i, min_int, max_int);
+	// test normal int <--> fp conversions
+	assert( fuzzy<double>::EQ( lm.fp_of_int(i), x) );
+	assert( lm.int_of_fp(x) == i );
+
+	// test fp --> int returning result as fp
+	if (i > min_int)
+		assert(
+		   fuzzy<double>::EQ(
+		      lm.fp_int_of_fp(x - 0.5*delta_fp), double(i) - 0.5
+				    )
+		      );
+	if (i < max_int)
+		assert(
+		   fuzzy<double>::EQ(
+		      lm.fp_int_of_fp(x + 0.5*delta_fp), double(i) + 0.5
+				    )
+		      );
+
+	// test silent rounding
+	if (i > min_int)
+		assert( lm.int_of_fp(x - 0.49*delta_fp, lm.round) == i );
+	if (i < max_int)
+		assert( lm.int_of_fp(x + 0.49*delta_fp, lm.round) == i );
+
+	// test floor/ceiling
+	if (i > min_int)
+		assert( lm.int_of_fp(x - 0.01*delta_fp, lm.floor) == i-1 );
+	if (i < max_int)
+		assert( lm.int_of_fp(x + 0.01*delta_fp, lm.ceiling) == i+1 );
+
+	// test zero-origin conversions
+	assert( lm.zero_origin_int(i) == i - min_int );
+	assert( lm.map_int(lm.zero_origin_int(i)) == i );
+	}
+
+printf("this should produce a pair of warning messages...\n");
+assert( lm.int_of_fp(min_fp + 0.01*delta_fp,
+		     linear_map<double>::warning) == min_int );
+assert( lm.int_of_fp(max_fp - 0.01*delta_fp,
+		     linear_map<double>::warning) == max_int );
+
+printf("this should produce an error message...\n");
+#ifdef TEST_RANGE
+assert( lm.int_of_fp(min_fp - 0.01*delta_fp) == min_int );
+#else
+assert( lm.int_of_fp(min_fp + 0.01*delta_fp) == min_int );
+#endif
+
+return 0;
+}