Index: trunk/tests/test_truncatedInverseSqrt.cc
===================================================================
--- trunk/tests/test_truncatedInverseSqrt.cc	(revision 0)
+++ trunk/tests/test_truncatedInverseSqrt.cc	(revision 571)
@@ -0,0 +1,50 @@
+#include "UnitTest++.h"
+
+#include "npstat/stat/InMemoryNtuple.hh"
+#include "npstat/stat/MultivariateSumsqAccumulator.hh"
+
+#include "npstat/rng/MersenneTwister.hh"
+
+#include "npstat/nm/truncatedInverseSqrt.hh"
+
+using namespace npstat;
+
+namespace {
+    TEST(truncatedInverseSqrt)
+    {
+        const double eps = 1.0e-12;
+
+        const unsigned dim = 5;
+        const unsigned npoints = 100;
+        const unsigned ntries = 3;
+
+        std::vector<double> row(dim);
+        MersenneTwister rng;
+
+        for (unsigned itry=0; itry<ntries; ++itry)
+        {
+            InMemoryNtuple<double> nt(simpleColumnNames(dim));
+
+            for (unsigned ipt=0; ipt<npoints; ++ipt)
+            {
+                for (unsigned i=0; i<dim; ++i)
+                    row[i] = rng();
+                nt.fill(&row[0], dim);
+            }
+
+            MultivariateSumAccumulator<> sums;
+            nt.cycleOverRows(sums);
+            MultivariateSumsqAccumulator<> sumsqs(sums);
+            nt.cycleOverRows(sumsqs);
+            const Matrix<double> covmat(sumsqs.covMat());
+
+            Matrix<double> invSqr;
+            const double frac = truncatedInverseSqrt(covmat, dim, &invSqr);
+            CHECK(frac == 0.0);
+
+            const Matrix<double>& invmat = covmat.symPDInv();
+            const double delta = (invSqr.TtimesThis() - invmat).maxAbsValue();
+            CHECK(delta < eps);
+        }
+    }
+}