initial import into cvs

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

8 files changed, 732 insertions, 0 deletions

diff --git a/src/jtutil/Makefile b/src/jtutil/Makefile
new file mode 100644
index 0000000..32f2384
--- /dev/null
+++ b/src/jtutil/Makefile
@@ -0,0 +1,95 @@
+# Makefile for my utility library (g++ + GNU make version)
+# $Id: Makefile,v 1.1 2001-06-15 11:21:34 jthorn Exp $
+#
+# Arguments:
+#	JT_GXX_FLAGS = (env in) My standard flags for g++(1).
+#	JT_LINTXX_FLAGS = (env in) My standard flags for lint++(1).
+#	LIB = (env in) My personal library directory.
+#	LIMBO = (env in) My personal "limbo" directory (for saving
+#			 old copies of libraries in).
+#
+# Targets:
+#	lib		==> lib-g lib-O
+#	lib-g		==> (re)build the library compiled with -g
+#			    and symlink the "public" library to the -g one
+#	lib-O		==> (re)build the library compiled with -O9
+#			    and symlink the "public" library to the -O one
+#	test		==> build test programs with -g
+#	lint		==> lint++ all the C files
+#	clean		==> delete object files, test drivers
+#	superclean	==> delete object files, test drivers, !library!
+#
+# Bugs:
+# - Dependencies on *.hh are omitted.
+# - All *.cc files (except test*) are assumed to go into the library.
+#
+
+###############################################################################
+
+CXX		:= g++
+LINTXX		:= lint++
+
+# default -- often overridden from command-line args
+CXXFLAGS :=
+
+CXXFLAGS-g	:= $(JT_GXX_FLAGS) $(CXXFLAGS) -g
+CXXFLAGS-O	:= $(JT_GXX_FLAGS) $(CXXFLAGS) -g -O9 -DNDEBUG
+LINTXXFLAGS	:= $(JT_LINTXX_FLAGS)
+
+LIBDIR		:= $(LIB)
+LIBUTILXX	:= libutil++.a
+LIBUTILXX-g	:= libutil++-g.a
+LIBUTILXX-O	:= libutil++-O.a
+
+# ... note '%' = wildcard character for filter() and filter-out()
+LIB_CCFILES	:= $(filter-out test%, $(wildcard *.cc))
+TEST_CCFILES	:= $(wildcard test*.cc)
+TEST_BINARIES	:= $(TEST_CCFILES:.cc=)
+
+###############################################################################
+
+.PHONY		: lib
+lib		: lib-g lib-O
+
+.PHONY		: lib-g
+lib-g		: $(LIBDIR)/$(LIBUTILXX-g)
+$(LIBDIR)/$(LIBUTILXX-g)	:: $(LIB_CCFILES)
+	$(CXX) -c $(CXXFLAGS-g) $?
+	ar rv $@ $(?:.cc=.o)
+	ranlib $@
+	rm $(?:.cc=.o)
+$(LIBDIR)/$(LIBUTILXX-g)	::
+	cd $(LIBDIR) && rm -f $(LIBUTILXX) && ln -s $(LIBUTILXX-g) $(LIBUTILXX)
+
+.PHONY		: lib-O
+lib-O		: $(LIBDIR)/$(LIBUTILXX-O)
+$(LIBDIR)/$(LIBUTILXX-O)	:: $(LIB_CCFILES)
+	$(CXX) -c $(CXXFLAGS-O) $?
+	ar rv $@ $(?:.cc=.o)
+	ranlib $@
+	rm $(?:.cc=.o)
+$(LIBDIR)/$(LIBUTILXX-O)	::
+	cd $(LIBDIR) && rm -f $(LIBUTILXX) && ln -s $(LIBUTILXX-O) $(LIBUTILXX)
+
+###############################################################################
+
+.PHONY		: test
+test		: $(TEST_BINARIES)
+$(TEST_BINARIES): %: %.cc
+	$(CXX) $(CXXFLAGS-g) -o $@ $? -lutil++-g -lutil-g -lm
+
+.PHONY		: lint
+lint		:
+	$(LINTXX) $(LINTXXFLAGS) $(LIB_CCFILES) $(TEST_CCFILES)
+
+.PHONY		: clean
+clean		:
+	-rm -f *.o $(TEST_BINARIES)
+
+.PHONY		: superclean
+superclean	: clean
+	-rm -f $(LIBDIR)/$(LIBUTILXX)
+	-test -f $(LIBDIR)/$(LIBUTILXX-g) \
+           && mv $(LIBDIR)/$(LIBUTILXX-g) $(LIMBO)
+	-test -f $(LIBDIR)/$(LIBUTILXX-O) \
+           && mv $(LIBDIR)/$(LIBUTILXX-O) $(LIMBO)
diff --git a/src/jtutil/README b/src/jtutil/README
new file mode 100644
index 0000000..d66b86c
--- /dev/null
+++ b/src/jtutil/README
@@ -0,0 +1,3 @@
+This directory holds low-level utility code that we use, but that's
+not really specific to this project -- things like min/max templates,
+1/2/3-dimensional array classes, fuzzy arithmetic routines, etc.
diff --git a/src/jtutil/error_exit.c b/src/jtutil/error_exit.c
new file mode 100644
index 0000000..26acdf7
--- /dev/null
+++ b/src/jtutil/error_exit.c
@@ -0,0 +1,48 @@
+/* error_exit -- print an error message and exit */
+/* $Id$ */
+/*
+ * error_exit -- print an error message and exit
+ */
+
+#include <stdio.h>
+#include <stdarg.h>
+#include <stdlib.h>
+
+#include "stdc.h"
+
+/*****************************************************************************/
+
+/*
+ * This function prints an error message (formatted via vfprintf(3S)),
+ * then exits.  It does *not* return to its caller.  Normally the exit
+ * is done via exit(2), but optionally it may be done via abort(2) to
+ * produce a core dump.  Despite not actually returning, this function
+ * is declared as returning an  int , so it may be easily used in
+ * conditional expressions like
+ *      foo = bar_ok ? baz(bar) : error_exit(...);
+ *
+ * Usage:
+ *	error_exit(exit_code, message, args...)
+ *
+ * ***** THIS VERSION IS FOR ANSI/ISO C *****
+ *
+ * Arguments:
+ * exit_code = (in) Exit code for exit(2), or PANIC_EXIT (defined in
+ *		    <jt/stdc.h>) to abort(2) instead.
+ * format = (in) vprintf(3S) format string for error message to be printed.
+ * args... = (in) Any additional arguments are (presumably) used in formatting
+ *		  the error message string.
+ */
+/*VARARGS*/
+int error_exit(const int exit_code, const char *const format, ...)
+{
+va_list ap;
+
+va_start(ap, format);
+vfprintf(stderr, format, ap);
+va_end(ap);
+
+if (exit_code == PANIC_EXIT)
+   then abort();						/*NOTREACHED*/
+   else exit(exit_code);					/*NOTREACHED*/
+}
diff --git a/src/jtutil/fuzzy.cc b/src/jtutil/fuzzy.cc
new file mode 100644
index 0000000..9706852
--- /dev/null
+++ b/src/jtutil/fuzzy.cc
@@ -0,0 +1,66 @@
+// fuzzy.cc -- template for fuzzy comparisons et al on floating point values
+// $Id$
+//
+// fuzzy::tolerance - comparison tolerance
+// fuzzy::EQ
+// fuzzy::is_integer
+// fuzzy::floor
+// fuzzy::ceiling
+//
+// *** instantiations of fuzzy template for <float> and <double>
+//
+
+#include "stdc.h"
+#include "util++.hh"
+
+//*****************************************************************************
+
+template <class fpt>
+bool fuzzy<fpt>::EQ(fpt x, fpt y)
+{
+fpt max_abs = jtutil::max(jtutil::abs(x), jtutil::abs(y));
+fpt epsilon = jtutil::max(tolerance, tolerance*max_abs);
+
+return jtutil::abs<fpt>(x-y) <= epsilon;
+}
+
+//*****************************************************************************
+
+template <class fpt>
+bool fuzzy<fpt>::is_integer(fpt x)
+{
+int i = round<fpt>::to_integer(x);
+return EQ(x, fpt(i));
+}
+
+//*****************************************************************************
+
+template <class fpt>
+int fuzzy<fpt>::floor(fpt x)
+{
+return fuzzy<fpt>::is_integer(x)
+       ? round<fpt>::to_integer(x)
+       : round<fpt>::floor(x);
+}
+
+//*****************************************************************************
+
+template <class fpt>
+int fuzzy<fpt>::ceiling(fpt x)
+{
+return fuzzy<fpt>::is_integer(x)
+       ? round<fpt>::to_integer(x)
+       : round<fpt>::ceiling(x);
+}
+
+//*****************************************************************************
+// *** instantiations of fuzzy template for <float> and <double>
+//*****************************************************************************
+
+// instantiation for <float>
+template class fuzzy<float>;
+float fuzzy<float>::tolerance = 1.0e-5;		// about 100 * FLT_EPSILON
+
+// instantiations for <double>
+template class fuzzy<double>;
+double fuzzy<double>::tolerance = 1.0e-12;	// about 1e4 * DBL_EPSILON
diff --git a/src/jtutil/round.cc b/src/jtutil/round.cc
new file mode 100644
index 0000000..3b1127f
--- /dev/null
+++ b/src/jtutil/round.cc
@@ -0,0 +1,66 @@
+// round.hh -- template for rounding floating point values
+// $Id$
+//
+// *** Implementation Notes ***
+// round::to_integer
+// round::floor
+// round::ceiling
+//
+// *** instantiations of round template for <float> and <double>
+//
+
+#include "stdc.h"
+#include "util++.hh"
+
+//*****************************************************************************
+// round.hh -- template for rounding floating point values
+//*****************************************************************************
+
+//
+// *** Implementation Notes ***
+//
+// We assume throughout this code that C++'s "built-in" conversion
+// from <fpt> to integer takes the floor, at least for zero or positive
+// values.
+//
+
+//*****************************************************************************
+
+// round to nearest integer, up for exact tie
+template <class fpt>
+int round<fpt>::to_integer(fpt x)
+{
+return (x >= 0.0)
+       ? int(x + 0.5)		// eg 3.6 --> int(4.1) = 4
+       : - int( (-x) + 0.5 );	// eg -3.6 --> - int(4.1) = -4
+}
+
+//*****************************************************************************
+
+template <class fpt>
+int round<fpt>::floor(fpt x)
+{
+return (x >= 0.0)
+       ? int(x)
+       : - ceiling(-x);
+}
+
+//*****************************************************************************
+
+template <class fpt>
+int round<fpt>::ceiling(fpt x)
+{
+return (x >= 0.0)
+       ? int(x) + (x != fpt(int(x)))
+       : - floor(-x);
+}
+
+//*****************************************************************************
+// *** instantiations of round template for <float> and <double>
+//*****************************************************************************
+
+// instantiation for <float>
+template class round<float>;
+
+// instantiations for <double>
+template class round<double>;
diff --git a/src/jtutil/test_fuzzy.cc b/src/jtutil/test_fuzzy.cc
new file mode 100644
index 0000000..68cd619
--- /dev/null
+++ b/src/jtutil/test_fuzzy.cc
@@ -0,0 +1,231 @@
+// test_fuzzy.cc -- test driver for fuzzy<float> and fuzzy<double>
+// $Id$
+//
+// main
+// check<fpt>
+// check_binary<fpt>
+// check_unary<fpt>
+// *** template instantiations
+
+#include <stdio.h>
+#include <assert.h>
+#include "stdc.h"
+#include "util++.hh"
+
+// prototypes
+template <class fpt>
+  void check(const char *print_string,
+	     fpt x, fpt y,
+	     char xy_relation,
+	     bool unary_flag,
+	     bool x_is_fuzzy_integer,
+	     int fuzzy_floor_of_x,
+	     int fuzzy_ceiling_of_x);
+template <class fpt>
+  void check_binary(fpt x, fpt y,
+		    char xy_relation);
+template <class fpt>
+  void check_unary(fpt x,
+		   bool x_is_fuzzy_integer,
+		   int fuzzy_floor_of_x,
+		   int fuzzy_ceiling_of_x);
+
+//*****************************************************************************
+
+int main(int argc, const char *const argv[])
+{
+const bool print_flag = (argc > 1);
+const char *double_str = print_flag ? "<double>" : NULL;
+const char *float_str  = print_flag ? "<float>" : NULL;
+
+	for (int ii = -10 ; ii <= +10 ; ++ii)
+	{
+	bool unary_flag = ((ii % 2) == 0);	// the "correct answers"
+						// for the unary tests
+						// assume f and d are integral
+
+	float f = 0.5*float(ii);
+	double d = 0.5*double(ii);
+	int i = ii / 2;
+
+	printf("testing %+-4g +/- various things...\n", d);
+
+
+	//                        x        y  relation  do unary  is_integer  floor  ceiling
+	//                                     (x,y)     tests?      (x)       (x)     (x)
+	check<double>(double_str, d - 0.49 , d,   '<',   unary_flag,  false,     i-1,    i  );
+	check<double>(double_str, d - 0.01 , d,   '<',   unary_flag,  false,     i-1,    i  );
+	check<double>(double_str, d - 1e-6 , d,   '<',   unary_flag,  false,     i-1,    i  );
+	check<double>(double_str, d - 1e-10, d,   '<',   unary_flag,  false,     i-1,    i  );
+	// ... distinct floating point numbers, but within tolerance
+	assert( d - 1e-14 != d );
+	check<double>(double_str, d - 1e-14, d,   '=',   unary_flag,  true ,     i  ,    i  );
+	check<double>(double_str, d        , d,   '=',   unary_flag,  true ,     i  ,    i  );
+	// ... distinct floating point numbers, but within tolerance
+	assert( d + 1e-14 != d );
+	check<double>(double_str, d + 1e-14, d,   '=',   unary_flag,  true ,     i  ,    i  );
+	check<double>(double_str, d + 1e-10, d,   '>',   unary_flag,  false,     i  ,    i+1);
+	check<double>(double_str, d + 1e-6 , d,   '>',   unary_flag,  false,     i  ,    i+1);
+	check<double>(double_str, d + 0.01 , d,   '>',   unary_flag,  false,     i  ,    i+1);
+	check<double>(double_str, d + 0.49 , d,   '>',   unary_flag,  false,     i  ,    i+1);
+
+	//                        x        y  relation  do unary  is_integer  floor  ceiling
+	//                                     (x,y)     tests?      (x)       (x)     (x)
+	check<float> (float_str , f - 0.49 , f,   '<',   unary_flag,  false,     i-1,    i  );
+	check<float> (float_str , f - 0.01 , f,   '<',   unary_flag,  false,     i-1,    i  );
+	check<float> (float_str , f - 1e-4 , f,   '<',   unary_flag,  false,     i-1,    i  );
+	// ... distinct floating point numbers, but within tolerance
+	assert( f - 1e-6 != f );
+	check<float> (float_str , f - 1e-6 , f,   '=',   unary_flag,  true ,     i  ,    i  );
+	// ... not enough precision to see as noninteger
+	check<float> (float_str , f - 1e-10, f,   '=',   unary_flag,  true ,     i  ,    i  );
+	check<float> (float_str , f        , f,   '=',   unary_flag,  true ,     i  ,    i  );
+	// ... not enough precision to see as noninteger
+	check<float> (float_str , f + 1e-10, f,   '=',   unary_flag,  true ,     i  ,    i  );
+	// ... distinct floating point numbers, but within tolerance
+	assert( f + 1e-6 != f );
+	check<float> (float_str , f + 1e-6 , f,   '=',   unary_flag,  true ,     i  ,    i  );
+	check<float> (float_str , f + 0.01 , f,   '>',   unary_flag,  false,     i  ,    i+1);
+	check<float> (float_str , f + 1e-4 , f,   '>',   unary_flag,  false,     i  ,    i+1);
+	check<float> (float_str , f + 0.49 , f,   '>',   unary_flag,  false,     i  ,    i+1);
+
+	if (print_flag)
+	   then printf("\n");
+	}
+
+printf("everything looks ok!\n");
+return 0;
+}
+
+//*****************************************************************************
+
+//
+// This function template is a driver to do all the tests for a single
+// argument or pair of arguments.  It calls  check_binary<fpt>()  to do the
+// binary-function tests, and optionally calls  check_unary<fpt>()  to do the
+// unary-function tests.
+//
+template <class fpt>
+  void check(const char *print_string,
+	     fpt x, fpt y,
+	     char xy_relation,
+	     bool unary_flag,
+	     bool x_is_fuzzy_integer,
+	     int fuzzy_floor_of_x,
+	     int fuzzy_ceiling_of_x)
+{
+if (print_string != NULL)
+   then printf("   testing %s %s: x=%.15f y=%.15f\n",
+	       print_string,
+	       unary_flag ? "op(x,y) + op(x)" : "op(x,y) only   ",
+	       double(x), double(y));
+
+check_binary<fpt>(x, y, xy_relation);
+
+if (unary_flag)
+   then check_unary(x, x_is_fuzzy_integer, fuzzy_floor_of_x, fuzzy_ceiling_of_x);
+}
+
+//*****************************************************************************
+
+//
+// This function template tests the binary relations
+//	fuzzy<fpt>::EQ()
+//	fuzzy<fpt>::NE()
+//	fuzzy<fpt>::LT()
+//	fuzzy<fpt>::GT()
+//	fuzzy<fpt>::LE()
+//	fuzzy<fpt>::GE()
+// for a single argument.
+//
+// The  xy_relation  argument must be one of '<', '=', or '>', indicating
+// the desired fuzzy relationship between the two arguments.  The function
+// infers the correct values of the binary functions from this.
+//
+template <class fpt>
+  void check_binary(fpt x, fpt y,
+		    char xy_relation)
+{
+//
+// basic properties which should always hold
+//
+
+// trichotomy (we have exactly one of <, =, >
+assert( fuzzy<fpt>::LT(x,y) + fuzzy<fpt>::EQ(x,y) + fuzzy<fpt>::GT(x,y) == 1 );
+
+// consistency of EQ() with NE()
+assert( fuzzy<fpt>::NE(x,y) == ! fuzzy<fpt>::EQ(x,y) );
+
+// consistency of LE() with LT() and EQ(), GE() with GT() and EQ()
+assert( fuzzy<fpt>::LE(x,y) == (fuzzy<fpt>::LT(x,y) || fuzzy<fpt>::EQ(x,y)) );
+assert( fuzzy<fpt>::GE(x,y) == (fuzzy<fpt>::GT(x,y) || fuzzy<fpt>::EQ(x,y)) );
+
+
+//
+// now check desired ordering relation
+// ... by virtue of the previous tests, we need only check LT(), EQ(), and GT()
+// ... in fact, by virtue of trichotomy (already tested), we really
+//     only need to assert() the true one for each case, but the slight
+//     redundancy of checking all of them seems nicer in this context
+//
+switch	(xy_relation)
+	{
+case '<':
+	assert( fuzzy<fpt>::LT(x,y) == true  );
+	assert( fuzzy<fpt>::EQ(x,y) == false );
+	assert( fuzzy<fpt>::GT(x,y) == false );
+	break;
+case '=':
+	assert( fuzzy<fpt>::LT(x,y) == false );
+	assert( fuzzy<fpt>::EQ(x,y) == true  );
+	assert( fuzzy<fpt>::GT(x,y) == false );
+	break;
+case '>':
+	assert( fuzzy<fpt>::LT(x,y) == false );
+	assert( fuzzy<fpt>::EQ(x,y) == false );
+	assert( fuzzy<fpt>::GT(x,y) == true  );
+	break;
+default:
+	jtuil::error_exit(PANIC_EXIT,
+"***** check_binary<fpt>: bad xy_relation=(int)'%c'\n"
+"                         (this should never happen!)\n"
+,
+			  int(xy_relation));				/*NOTREACHED*/
+	}
+}
+
+//*****************************************************************************
+
+//
+// This function template tests the unary functions
+//	fuzzy<fpt>::is_integer()
+//	fuzzy<fpt>::floor()
+//	fuzzy<fpt>::ceiling()
+// for a single argument.
+//
+template <class fpt>
+  void check_unary(fpt x,
+		   // remaining arguments are the correct results of...
+		   bool x_is_fuzzy_integer,	// ... fuzzy<fpt>::is_integer(x)
+		   int fuzzy_floor_of_x,	// ... fuzzy<fpt>::floor(x)
+		   int fuzzy_ceiling_of_x)	// ... fuzzy<fpt>::ceiling(x)
+{
+assert( fuzzy<fpt>::is_integer(x) == x_is_fuzzy_integer );
+assert( fuzzy<fpt>::floor(x) == fuzzy_floor_of_x );
+assert( fuzzy<fpt>::ceiling(x) == fuzzy_ceiling_of_x );
+}
+
+//*****************************************************************************
+
+//
+// template instantiations for <float> and <double>
+//
+
+template void check<float>(const char *, float, float, char, bool, bool, int, int);
+template void check<double>(const char *, double, double, char, bool, bool, int, int);
+
+template void check_binary<float>(float, float, char);
+template void check_binary<double>(double, double, char);
+
+template void check_unary<float>(float, bool, int, int);
+template void check_unary<double>(double, bool, int, int);
diff --git a/src/jtutil/test_round.cc b/src/jtutil/test_round.cc
new file mode 100644
index 0000000..075209b
--- /dev/null
+++ b/src/jtutil/test_round.cc
@@ -0,0 +1,104 @@
+// test_round.cc -- test driver for round<fpt> and round<double>
+// $Id$
+
+#include <stdio.h>
+#include <assert.h>
+#include "stdc.h"
+#include "util++.hh"
+
+//
+// Usage:
+//	test_round x
+//
+int main(int argc, const char *const argv[])
+{
+	for (int i = -10 ; i <= +10 ; ++i)
+	{
+	printf("testing %d +/- various things...\n", i);
+
+	float f = float(i);
+	double d = double(i);
+
+
+	//
+	// round<fpt>::to_integer()
+	//
+
+	assert( round<double>::to_integer(d - 0.49) == i );
+	assert( round<double>::to_integer(d - 0.01) == i );
+	assert( round<double>::to_integer(d - 1e-4) == i );
+	assert( round<double>::to_integer(d - 1e-10) == i );
+	assert( round<double>::to_integer(d        ) == i );
+	assert( round<double>::to_integer(d + 1e-10) == i );
+	assert( round<double>::to_integer(d + 1e-4) == i );
+	assert( round<double>::to_integer(d + 0.01) == i );
+	assert( round<double>::to_integer(d + 0.49) == i );
+
+	assert( round<float>::to_integer(f - 0.49) == i );
+	assert( round<float>::to_integer(f - 0.01) == i );
+	assert( round<float>::to_integer(f - 1e-4) == i );
+	assert( round<float>::to_integer(f - 1e-10) == i );
+	assert( round<float>::to_integer(f        ) == i );
+	assert( round<float>::to_integer(f + 1e-10) == i );
+	assert( round<float>::to_integer(f + 1e-4) == i );
+	assert( round<float>::to_integer(f + 0.01) == i );
+	assert( round<float>::to_integer(f + 0.49) == i );
+
+
+	//
+	// round<fpt>::floor()
+	//
+
+	assert( round<double>::floor(d - 0.49) == i-1 );
+	assert( round<double>::floor(d - 0.01) == i-1);
+	assert( round<double>::floor(d - 1e-4) == i-1 );
+	assert( round<double>::floor(d - 1e-10) == i-1 );
+	assert( round<double>::floor(d        ) == i );
+	assert( round<double>::floor(d + 1e-10) == i );
+	assert( round<double>::floor(d + 1e-4) == i );
+	assert( round<double>::floor(d + 0.01) == i );
+	assert( round<double>::floor(d + 0.49) == i );
+
+	assert( round<float>::floor(f - 0.49) == i-1 );
+	assert( round<float>::floor(f - 0.01) == i-1 );
+	assert( round<float>::floor(f - 1e-4) == i-1 );
+	assert( round<float>::floor(f - 1e-10) == ((i == 0) ? i-1 : i) );
+					// i != 0 ==> not enough precision
+					//            to see as noninteger
+	assert( round<float>::floor(f        ) == i );
+	assert( round<float>::floor(f + 1e-10) == i );
+	assert( round<float>::floor(f + 1e-4) == i );
+	assert( round<float>::floor(f + 0.01) == i );
+	assert( round<float>::floor(f + 0.49) == i );
+
+
+	//
+	// round<fpt>::ceiling()
+	//
+
+	assert( round<double>::ceiling(d - 0.49) == i );
+	assert( round<double>::ceiling(d - 0.01) == i);
+	assert( round<double>::ceiling(d - 1e-4) == i );
+	assert( round<double>::ceiling(d - 1e-10) == i );
+	assert( round<double>::ceiling(d        ) == i );
+	assert( round<double>::ceiling(d + 1e-10) == i+1 );
+	assert( round<double>::ceiling(d + 1e-4) == i+1 );
+	assert( round<double>::ceiling(d + 0.01) == i+1 );
+	assert( round<double>::ceiling(d + 0.49) == i+1 );
+
+	assert( round<float>::ceiling(f - 0.49) == i );
+	assert( round<float>::ceiling(f - 0.01) == i);
+	assert( round<float>::ceiling(f - 1e-4) == i );
+	assert( round<float>::ceiling(f - 1e-10) == i );
+	assert( round<float>::ceiling(f        ) == i );
+	assert( round<float>::ceiling(f + 1e-10) == ((i == 0) ? i+1 : i) );
+					// i != 0 ==> not enough precision
+					//            to see as noninteger
+	assert( round<float>::ceiling(f + 1e-4) == i+1 );
+	assert( round<float>::ceiling(f + 0.01) == i+1 );
+	assert( round<float>::ceiling(f + 0.49) == i+1 );
+	}
+
+printf("everything looks ok!\n");
+return 0;
+}
diff --git a/src/jtutil/util++2.hh b/src/jtutil/util++2.hh
new file mode 100644
index 0000000..9c37159
--- /dev/null
+++ b/src/jtutil/util++2.hh
@@ -0,0 +1,119 @@
+// util++.hh -- stuff for my C++ utility library
+// $Id$
+//
+// jtutil::{min,max} - min/max templates
+// jtutil::abs - absolute value template
+// jtutil::pow* - raise floating point value to small-integer power
+// fuzzy::<fpt> - fuzzy floating point comparisons and other operations
+// round::<fpt> - floating point rounding
+//
+
+//******************************************************************************
+
+namespace jtutil
+	  {
+
+// how many integers are in the closed interval [low,high]
+inline int how_many_in_range(int low, int high)
+	{ return high - low + 1; }
+
+//
+// minimum/maximum templates (valid for both integer and floating-point types)
+// FIXME: these are supposed to be in the STL,
+//	  but in practice they often aren't (yet) :(
+//	  eg the gcc 2.95.2 library doesn't have them :(
+//
+template <typename num>
+  inline num min(num x, num y) { return x < y ? x : y; };
+template <typename num>
+  inline num max(num x, num y) { return x > y ? x : y; };
+
+//
+// absolute value template
+// FIXME: this ought to be somewhere in the STL?
+//
+template <typename fpt>
+  inline fpt abs(fpt x) { return (x > 0.0) ? x : -x; }
+
+//
+// These functions raise their arguments (presumably floating-point)
+// to various small-integer powers.
+// FIXME: do we ever use these?
+//
+template <typename fpt>
+  inline fpt pow2(fpt x) { return x*x; }
+template <typename fpt>
+  inline fpt pow3(fpt x) { return x*x*x; }
+template <typename fpt>
+  inline fpt pow4(fpt x) { return pow2(pow2(x)); }
+template <typename fpt>
+  inline fpt pow5(fpt x) { return x * pow4(x); }
+template <typename fpt>
+  inline fpt pow6(fpt x) { return pow3(pow2(x)); }
+template <typename fpt>
+  inline fpt pow7(fpt x) { return x * pow6(x); }
+template <typename fpt>
+  inline fpt pow8(fpt x) { return pow2(pow2(pow2(x))); }
+
+	  // close namespace jtutil::
+	  };
+
+//******************************************************************************
+
+//
+// This template does fuzzy comparisons and related operations on
+// floating point values, parameterized by the floating point type.
+//
+// The fuzzy comparison semantics are based on those of APL, but are
+// modified to use an absolute error tolerance for values close to 0.
+//
+
+// this template class has only static members
+// ... it's a class, not a namespace, because we want to express the
+//     semantics that the entire class is a single template, rather
+//     than the individual members being conceptually-unrelated templates
+// ... moreover, we need the *data* member (template)  tolerance , and
+//     it seems C++ doesn't grok data templates which aren't in a class
+template <class fpt>
+class	fuzzy
+	{
+public:
+	// comparison tolerance
+	// ... must be explicitly initialized when instantiating
+	//     for a new <fpt> type, see "fuzzy.cc" for details/examples
+	// ... may be modified by user code if needed
+	static fpt tolerance;
+
+	// fuzzy commparisons
+	static bool EQ(fpt x, fpt y);
+	static bool NE(fpt x, fpt y) { return ! EQ(x,y); }
+	static bool LT(fpt x, fpt y) { return EQ(x,y) ? false : (x < y); }
+	static bool LE(fpt x, fpt y) { return EQ(x,y) ? true  : (x < y); }
+	static bool GT(fpt x, fpt y) { return EQ(x,y) ? false : (x > y); }
+	static bool GE(fpt x, fpt y) { return EQ(x,y) ? true  : (x > y); }
+
+	static bool is_integer(fpt x);	// is x fuzzily an integer?
+	static int floor(fpt x);	// round x fuzzily down to integer
+	static int ceiling(fpt x);	// round x fuzzily up to integer
+	};
+
+//******************************************************************************
+
+//
+// This template does machine-independent rounding of floating point
+// values, parameterized by the floating point type.
+//
+
+// this template class has only static members
+// ... it's a class, not a namespace, because we want to express the
+//     semantics that the entire class is a single template, rather
+//     than the individual members being conceptually-unrelated templates
+template <class fpt>
+class	round
+	{
+public:
+	static int to_integer(fpt x);	// round to nearest integer
+
+	static int floor(fpt x);	// round down to integer
+	static int ceiling(fpt x);	// round up to integer
+	};