Mercurial > octave-nkf

changeset 8335:64cf956a109c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
templatize & fix DET
author Jaroslav Hajek <highegg@gmail.com>
date Wed, 19 Nov 2008 11:23:07 +0100
parents 70dd33450061
children 9813c07ca946
files liboctave/CMatrix.cc liboctave/CMatrix.h liboctave/CSparse.cc liboctave/CSparse.h liboctave/ChangeLog liboctave/CmplxDET.cc liboctave/CmplxDET.h liboctave/DET.h liboctave/Makefile.in liboctave/dMatrix.cc liboctave/dMatrix.h liboctave/dSparse.cc liboctave/dSparse.h liboctave/dbleDET.cc liboctave/dbleDET.h liboctave/fCMatrix.cc liboctave/fCMatrix.h liboctave/fCmplxDET.cc liboctave/fCmplxDET.h liboctave/fMatrix.cc liboctave/fMatrix.h liboctave/floatDET.cc liboctave/floatDET.h liboctave/mx-defs.h liboctave/mx-ext.h src/ChangeLog src/DLD-FUNCTIONS/det.cc
diffstat 27 files changed, 158 insertions(+), 964 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CMatrix.cc	Tue Nov 18 22:55:13 2008 +0100
+++ b/liboctave/CMatrix.cc	Wed Nov 19 11:23:07 2008 +0100
@@ -40,7 +40,7 @@
 #include "Array-util.h"
 #include "CMatrix.h"
 #include "CmplxAEPBAL.h"
-#include "CmplxDET.h"
+#include "DET.h"
 #include "CmplxSCHUR.h"
 #include "CmplxSVD.h"
 #include "CmplxCHOL.h"
@@ -1631,33 +1631,13 @@
 	    } 
 	  else 
 	    {
-	      Complex c = 1.0;
-	      int e = 0;
-
+              retval = ComplexDET (1.0);
+              
 	      for (octave_idx_type i = 0; i < nc; i++) 
 		{
-		  if (ipvt(i) != (i+1))
-		    c = -c;
-
-		  c *= atmp(i,i);
-
-		  if (c == 0.0)
-		    break;
-
-		  while (std::abs(c) < 0.5)
-		    {
-		      c *= 2.0;
-		      e--;
-		    }
-
-		  while (std::abs(c) >= 2.0)
-		    {
-		      c /= 2.0;
-		      e++;
-		    }
-		}
-
-	      retval = ComplexDET (c, e);
+                  Complex c = atmp(i,i);
+                  retval *= (ipvt(i) != (i+1)) ? -c : c;
+                }
 	    }
 	}
     }
--- a/liboctave/CMatrix.h	Tue Nov 18 22:55:13 2008 +0100
+++ b/liboctave/CMatrix.h	Wed Nov 19 11:23:07 2008 +0100
@@ -31,6 +31,7 @@
 #include "mx-defs.h"
 #include "mx-op-defs.h"
 #include "oct-cmplx.h"
+#include "DET.h"
 
 class
 OCTAVE_API
--- a/liboctave/CSparse.cc	Tue Nov 18 22:55:13 2008 +0100
+++ b/liboctave/CSparse.cc	Wed Nov 19 11:23:07 2008 +0100
@@ -1114,10 +1114,7 @@
 
   if (nr == 0 || nc == 0 || nr != nc)
     {
-      Complex d[2];
-      d[0] = 1.0;
-      d[1] = 0.0;
-      retval = ComplexDET (d);
+      retval = ComplexDET (1.0);
     }
   else
     {
@@ -1201,13 +1198,10 @@
 	    {
 	      UMFPACK_ZNAME (report_numeric) (Numeric, control);
 
-	      Complex d[2];
-	      double d_exponent;
-
-	      status = UMFPACK_ZNAME (get_determinant) 
-		(reinterpret_cast<double *> (&d[0]), 0, &d_exponent,
-		 Numeric, info);
-	      d[1] = d_exponent;
+	      double c10[2], e10;
+
+              status = UMFPACK_ZNAME (get_determinant) (c10, 0, &e10,
+                                                        Numeric, info);
 
 	      if (status < 0)
 		{
@@ -1218,7 +1212,7 @@
 		  UMFPACK_ZNAME (report_info) (control, info);
 		}
 	      else
-		retval = ComplexDET (d);
+		retval = ComplexDET (Complex (c10[0], c10[1]), e10, 10);
 		  
 	      UMFPACK_ZNAME (free_numeric) (&Numeric);
 	    }
--- a/liboctave/CSparse.h	Tue Nov 18 22:55:13 2008 +0100
+++ b/liboctave/CSparse.h	Wed Nov 19 11:23:07 2008 +0100
@@ -31,7 +31,7 @@
 #include "CColVector.h"
 #include "oct-cmplx.h"
 
-#include "CmplxDET.h"
+#include "DET.h"
 #include "MSparse.h"
 #include "MSparse-defs.h"
 #include "Sparse-op-defs.h"
--- a/liboctave/ChangeLog	Tue Nov 18 22:55:13 2008 +0100
+++ b/liboctave/ChangeLog	Wed Nov 19 11:23:07 2008 +0100
@@ -1,3 +1,19 @@
+2008-11-19  Jaroslav Hajek  <highegg@gmail.com>
+
+	* DET.h: New source.
+	* CmplxDET.cc, CmplxDET.h, dbleDET.cc, dbleDET.h, fCmplxDET.cc,
+	fCmplxDET.h, floatDET.cc, floatDET.h: Remove.
+	* Makefile.in: Reflect changes.
+	* mx-defs.h: Remove DET decls.
+	* mx-ext.h, dMatrix.h, fMatrix.h, CMatrix.h, fCMatrix.h,
+	dSparse.h, CSparse.h: Include only DET.h.
+	* dMatrix.cc (Matrix::determinant),
+	fMatrix.cc (FloatMatrix::determinant),
+	CMatrix.cc (ComplexMatrix::determinant),
+	fCMatrix.cc (FloatComplexMatrix::determinant),
+	dSparse.cc (SparseMatrix::determinant),
+	CSparse.cc (SparseComplexMatrix::determinant): Use new class.
+
 2008-11-18  David Bateman  <dbateman@free.fr>
 
 	* file-ops.cc (std::string file_ops::tilde_expand (const 
--- a/liboctave/CmplxDET.cc	Tue Nov 18 22:55:13 2008 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,86 +0,0 @@
-/*
-
-Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2006,
-              2007 John W. Eaton
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <cassert>
-#include <cfloat>
-
-#include "CmplxDET.h"
-#include "lo-mappers.h"
-#include "lo-math.h"
-#include "oct-cmplx.h"
-
-bool
-ComplexDET::value_will_overflow (void) const
-{
-  return base2
-    ? (e2 + 1 > xlog2 (DBL_MAX) ? 1 : 0)
-    : (e10 + 1 > log10 (DBL_MAX) ? 1 : 0);
-}
-
-bool
-ComplexDET::value_will_underflow (void) const
-{
-  return base2
-    ? (e2 - 1 < xlog2 (DBL_MIN) ? 1 : 0)
-    : (e10 - 1 < log10 (DBL_MIN) ? 1 : 0);
-}
-
-void
-ComplexDET::initialize10 (void)
-{
-  if (c2 != 0.0)
-    {
-      double etmp = e2 / xlog2 (static_cast<double>(10));
-      e10 = static_cast<int> (xround (etmp));
-      etmp -= e10;
-      c10 = c2 * pow (10.0, etmp);
-    }
-}
-
-void
-ComplexDET::initialize2 (void)
-{
-  if (c10 != 0.0)
-    {
-      double etmp = e10 / log10 (2.0);
-      e2 = static_cast<int> (xround (etmp));
-      etmp -= e2;
-      c2 = c10 * xexp2 (etmp);
-    }
-}
-
-Complex
-ComplexDET::value (void) const
-{
-  return base2 ? c2 * xexp2 (static_cast<double>(e2)) : c10 * pow (10.0, e10);
-}
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
--- a/liboctave/CmplxDET.h	Tue Nov 18 22:55:13 2008 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,121 +0,0 @@
-/*
-
-Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
-              2007 John W. Eaton
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#if !defined (octave_ComplexDET_h)
-#define octave_ComplexDET_h 1
-
-#include <iostream>
-
-#include "oct-cmplx.h"
-
-// FIXME -- we could use templates here; compare with dbleDET.h
-
-class
-OCTAVE_API
-ComplexDET
-{
-friend class ComplexMatrix;
-friend class SparseComplexMatrix;
-
-public:
-
-  ComplexDET (void) : c2 (0), c10 (0), e2 (0), e10 (0), base2 (false) { }
-
-  ComplexDET (const ComplexDET& a)
-    : c2 (a.c2), c10 (a.c10), e2 (a.e2), e10 (a.e10), base2 (a.base2)
-    { }
-
-  ComplexDET& operator = (const ComplexDET& a)
-    {
-      if (this != &a)
-	{
-	  c2 = a.c2;
-	  e2 = a.e2;
-
-	  c10 = a.c10;
-	  e10 = a.e10;
-
-	  base2 = a.base2;
-	}
-      return *this;
-    }
-
-  bool value_will_overflow (void) const;
-  bool value_will_underflow (void) const;
-
-  // These two functions were originally defined in base 10, so we are
-  // preserving that interface here.
-
-  Complex coefficient (void) const { return coefficient10 (); }
-  int exponent (void) const { return exponent10 (); }
-
-  Complex coefficient10 (void) const { return c10; }
-  int exponent10 (void) const { return e10; }
-
-  Complex coefficient2 (void) const { return c2; }
-  int exponent2 (void) const { return e2; }
-
-  Complex value (void) const;
-
-  friend std::ostream&  operator << (std::ostream& os, const ComplexDET& a);
-
-private:
-
-  // Constructed this way, we assume base 2.
-
-  ComplexDET (const Complex& c, int e)
-    : c2 (c), c10 (0), e2 (e), e10 (0), base2 (true)
-    {
-      initialize10 ();
-    }
-
-  // Original interface had only this constructor and it was assumed
-  // to be base 10, so we are preserving that interface here.
-
-  ComplexDET (const Complex *d)
-    : c2 (0), c10 (d[0]), e2 (0), e10 (static_cast<int> (d[1].real ())),
-      base2 (false)
-    {
-      initialize2 ();
-    }
-
-  void initialize2 (void);
-  void initialize10 (void);
-
-  Complex c2;
-  Complex c10;
-
-  int e2;
-  int e10;
-
-  // TRUE means the original values were provided in base 2.
-  bool base2;
-};
-
-#endif
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/DET.h	Wed Nov 19 11:23:07 2008 +0100
@@ -0,0 +1,87 @@
+/*
+
+Copyright (C) 2008 Jaroslav Hajek
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#if !defined (octave_DET_h)
+#define octave_DET_h 1
+
+#include <cmath>
+#include "oct-cmplx.h"
+#include "lo-mappers.h"
+
+template <class T>
+class
+OCTAVE_API
+base_det
+{
+public:
+
+  base_det (T c = 0, int e = 0) 
+    { 
+      c2 = xlog2 (c, e2); 
+      e2 += e; 
+    }
+
+  base_det (T c, double e, double b) 
+    { 
+      e *= xlog2 (b);
+      e2 = e;
+      c *= xexp2 (e - e2);
+      int f;
+      c2 = xlog2 (c, f);
+      e2 += f;
+    }
+
+  base_det (const base_det& a) : c2(a.c2), e2(a.e2) { }
+
+  base_det& operator = (const base_det& a)
+    {
+      c2 = a.c2;
+      e2 = a.e2;
+      return *this;
+    }
+
+  T coef (void) const { return c2; }
+  int exp (void) const { return e2; }
+
+  T value () const { return c2 * static_cast<T> (std::ldexp (1.0, e2)); }
+  operator T () const { return value (); }
+
+  void operator *= (T t)
+    {
+      int e;
+      c2 *= xlog2 (t, e);
+      e2 += e;
+    }
+
+private:
+
+  T c2;
+  int e2;
+};
+
+// Provide the old types by typedefs.
+typedef base_det<double> DET;
+typedef base_det<float> FloatDET;
+typedef base_det<Complex> ComplexDET;
+typedef base_det<FloatComplex> FloatComplexDET;
+
+#endif
--- a/liboctave/Makefile.in	Tue Nov 18 22:55:13 2008 +0100
+++ b/liboctave/Makefile.in	Wed Nov 19 11:23:07 2008 +0100
@@ -46,10 +46,10 @@
 	base-lu.h dim-vector.h mx-base.h mx-op-defs.h \
 	mx-defs.h mx-ext.h CColVector.h CDiagMatrix.h CMatrix.h \
 	CNDArray.h CRowVector.h CmplxAEPBAL.h CmplxCHOL.h \
-	CmplxDET.h CmplxGEPBAL.h CmplxHESS.h CmplxLU.h CmplxQR.h CmplxQRP.h \
+	CmplxGEPBAL.h CmplxHESS.h CmplxLU.h CmplxQR.h CmplxQRP.h \
 	CmplxSCHUR.h CmplxSVD.h EIG.h fEIG.h boolMatrix.h boolNDArray.h \
 	chMatrix.h chNDArray.h dColVector.h dDiagMatrix.h dMatrix.h \
-	dNDArray.h dRowVector.h dbleAEPBAL.h dbleCHOL.h dbleDET.h \
+	dNDArray.h dRowVector.h dbleAEPBAL.h dbleCHOL.h DET.h \
 	dbleGEPBAL.h dbleHESS.h dbleLU.h dbleQR.h dbleQRP.h dbleSCHUR.h \
 	dbleSVD.h boolSparse.h CSparse.h dSparse.h MSparse-defs.h MSparse.h \
 	Sparse.h sparse-base-lu.h SparseCmplxLU.h SparsedbleLU.h \
@@ -62,9 +62,9 @@
 	fCColVector.h fCRowVector.h fCDiagMatrix.h fCMatrix.h fCNDArray.h \
 	fColVector.h fRowVector.h fDiagMatrix.h fMatrix.h fNDArray.h \
 	fCmplxAEPBAL.h fCmplxGEPBAL.h fCmplxHESS.h fCmplxCHOL.h \
-	fCmplxDET.h fCmplxLU.h fCmplxSCHUR.h fCmplxSVD.h fCmplxQR.h \
+	fCmplxLU.h fCmplxSCHUR.h fCmplxSVD.h fCmplxQR.h \
 	fCmplxQRP.h floatAEPBAL.h \
-	floatCHOL.h floatDET.h floatGEPBAL.h floatHESS.h floatLU.h \
+	floatCHOL.h floatGEPBAL.h floatHESS.h floatLU.h \
 	floatSCHUR.h floatSVD.h floatQR.h floatQRP.h
 
 MX_OP_INC := $(shell $(AWK) -f $(srcdir)/mk-ops.awk prefix=mx list_h_files=1 $(srcdir)/mx-ops)
@@ -116,12 +116,12 @@
 
 MATRIX_SRC := Array-util.cc CColVector.cc \
 	CDiagMatrix.cc CMatrix.cc CNDArray.cc CRowVector.cc \
-	CmplxAEPBAL.cc CmplxCHOL.cc CmplxDET.cc CmplxGEPBAL.cc CmplxHESS.cc \
+	CmplxAEPBAL.cc CmplxCHOL.cc CmplxGEPBAL.cc CmplxHESS.cc \
 	CmplxLU.cc CmplxQR.cc CmplxQRP.cc CmplxSCHUR.cc CmplxSVD.cc \
 	EIG.cc fEIG.cc boolMatrix.cc boolNDArray.cc chMatrix.cc \
 	chNDArray.cc dColVector.cc dDiagMatrix.cc dMatrix.cc \
 	dNDArray.cc dRowVector.cc dbleAEPBAL.cc dbleCHOL.cc \
-	dbleDET.cc dbleGEPBAL.cc dbleHESS.cc dbleLU.cc dbleQR.cc dbleQRP.cc \
+	dbleGEPBAL.cc dbleHESS.cc dbleLU.cc dbleQR.cc dbleQRP.cc \
 	dbleSCHUR.cc dbleSVD.cc boolSparse.cc CSparse.cc dSparse.cc \
 	MSparse.cc Sparse.cc SparseCmplxLU.cc SparsedbleLU.cc \
 	SparseCmplxCHOL.cc SparsedbleCHOL.cc \
@@ -130,9 +130,9 @@
 	int32NDArray.cc uint32NDArray.cc int64NDArray.cc uint64NDArray.cc \
 	fCColVector.cc fCRowVector.cc fCDiagMatrix.cc fCMatrix.cc fCNDArray.cc \
 	fColVector.cc fRowVector.cc fDiagMatrix.cc fMatrix.cc fNDArray.cc \
-	fCmplxAEPBAL.cc fCmplxCHOL.cc fCmplxDET.cc fCmplxGEPBAL.cc \
+	fCmplxAEPBAL.cc fCmplxCHOL.cc fCmplxGEPBAL.cc \
 	fCmplxHESS.cc fCmplxLU.cc fCmplxSCHUR.cc fCmplxSVD.cc fCmplxQR.cc \
-	fCmplxQRP.cc floatAEPBAL.cc floatCHOL.cc floatDET.cc \
+	fCmplxQRP.cc floatAEPBAL.cc floatCHOL.cc \
 	floatGEPBAL.cc floatHESS.cc floatLU.cc \
 	floatSCHUR.cc floatSVD.cc floatQR.cc floatQRP.cc
 
--- a/liboctave/dMatrix.cc	Tue Nov 18 22:55:13 2008 +0100
+++ b/liboctave/dMatrix.cc	Wed Nov 19 11:23:07 2008 +0100
@@ -36,7 +36,7 @@
 #include "byte-swap.h"
 #include "dMatrix.h"
 #include "dbleAEPBAL.h"
-#include "dbleDET.h"
+#include "DET.h"
 #include "dbleSCHUR.h"
 #include "dbleSVD.h"
 #include "dbleCHOL.h"
@@ -1297,33 +1297,13 @@
 	    } 
 	  else 
 	    {
-	      double c = 1.0;
-	      int e = 0;
-
+              retval = DET (1.0);
+              
 	      for (octave_idx_type i = 0; i < nc; i++) 
 		{
-		  if (ipvt(i) != (i+1))
-		    c = -c;
-
-		  c *= atmp(i,i);
-
-		  if (c == 0.0)
-		    break;
-
-		  while (fabs (c) < 0.5)
-		    {
-		      c *= 2.0;
-		      e--;
-		    }
-
-		  while (fabs (c) >= 2.0)
-		    {
-		      c /= 2.0;
-		      e++;
-		    }
-		}
-
-	      retval = DET (c, e);
+                  double c = atmp(i,i);
+                  retval *= (ipvt(i) != (i+1)) ? -c : c;
+                }
 	    }
 	}
     }
--- a/liboctave/dMatrix.h	Tue Nov 18 22:55:13 2008 +0100
+++ b/liboctave/dMatrix.h	Wed Nov 19 11:23:07 2008 +0100
@@ -30,6 +30,7 @@
 
 #include "mx-defs.h"
 #include "mx-op-defs.h"
+#include "DET.h"
 
 class
 OCTAVE_API
--- a/liboctave/dSparse.cc	Tue Nov 18 22:55:13 2008 +0100
+++ b/liboctave/dSparse.cc	Wed Nov 19 11:23:07 2008 +0100
@@ -1188,10 +1188,7 @@
 
   if (nr == 0 || nc == 0 || nr != nc)
     {
-      double d[2];
-      d[0] = 1.0;
-      d[1] = 0.0;
-      retval = DET (d);
+      retval = DET (1.0);
     }
   else
     {
@@ -1270,10 +1267,9 @@
 	    {
 	      UMFPACK_DNAME (report_numeric) (Numeric, control);
 
-	      double d[2];
-
-	      status = UMFPACK_DNAME (get_determinant) (&d[0],
-						   &d[1], Numeric, info);
+	      double c10, e10;
+
+	      status = UMFPACK_DNAME (get_determinant) (&c10, &e10, Numeric, info);
 
 	      if (status < 0)
 		{
@@ -1284,7 +1280,7 @@
 		  UMFPACK_DNAME (report_info) (control, info);
 		}
 	      else
-		retval = DET (d);
+		retval = DET (c10, e10, 10);
 
 	      UMFPACK_DNAME (free_numeric) (&Numeric);
 	    }
--- a/liboctave/dSparse.h	Tue Nov 18 22:55:13 2008 +0100
+++ b/liboctave/dSparse.h	Wed Nov 19 11:23:07 2008 +0100
@@ -30,7 +30,7 @@
 #include "dColVector.h"
 #include "CColVector.h"
 
-#include "dbleDET.h"
+#include "DET.h"
 #include "MSparse.h"
 #include "MSparse-defs.h"
 #include "Sparse-op-defs.h"
--- a/liboctave/dbleDET.cc	Tue Nov 18 22:55:13 2008 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,84 +0,0 @@
-/*
-
-Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2006,
-              2007 John W. Eaton
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <cfloat>
-
-#include "dbleDET.h"
-#include "lo-mappers.h"
-#include "lo-math.h"
-
-bool
-DET::value_will_overflow (void) const
-{
-  return base2
-    ? (e2 + 1 > xlog2 (DBL_MAX) ? 1 : 0)
-    : (e10 + 1 > log10 (DBL_MAX) ? 1 : 0);
-}
-
-bool
-DET::value_will_underflow (void) const
-{
-  return base2
-    ? (e2 - 1 < xlog2 (DBL_MIN) ? 1 : 0)
-    : (e10 - 1 < log10 (DBL_MIN) ? 1 : 0);
-}
-
-void
-DET::initialize10 (void)
-{
-  if (c2 != 0.0)
-    {
-      double etmp = e2 / xlog2 (static_cast<double>(10));
-      e10 = static_cast<int> (xround (etmp));
-      etmp -= e10;
-      c10 = c2 * pow (10.0, etmp);
-    }
-}
-
-void
-DET::initialize2 (void)
-{
-  if (c10 != 0.0)
-    {
-      double etmp = e10 / log10 (2.0);
-      e2 = static_cast<int> (xround (etmp));
-      etmp -= e2;
-      c2 = c10 * xexp2 (etmp);
-    }
-}
-
-double
-DET::value (void) const
-{
-  return base2 ? c2 * xexp2 (static_cast<double>(e2)) : c10 * pow (10.0, e10);
-}
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
--- a/liboctave/dbleDET.h	Tue Nov 18 22:55:13 2008 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,118 +0,0 @@
-/*
-
-Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
-              2007 John W. Eaton
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#if !defined (octave_DET_h)
-#define octave_DET_h 1
-
-#include <iostream>
-
-// FIXME -- we could use templates here; compare with CmplxDET.h
-
-class
-OCTAVE_API
-DET
-{
-friend class Matrix;
-friend class SparseMatrix;
-
-public:
-
-  DET (void) : c2 (0), c10 (0), e2 (0), e10 (0), base2 (false) { }
-
-  DET (const DET& a)
-    : c2 (a.c2), c10 (a.c10), e2 (a.e2), e10 (a.e10), base2 (a.base2)
-    { }
-
-  DET& operator = (const DET& a)
-    {
-      if (this != &a)
-	{
-	  c2 = a.c2;
-	  e2 = a.e2;
-
-	  c10 = a.c10;
-	  e10 = a.e10;
-
-	  base2 = a.base2;
-	}
-      return *this;
-    }
-
-  bool value_will_overflow (void) const;
-  bool value_will_underflow (void) const;
-
-  // These two functions were originally defined in base 10, so we are
-  // preserving that interface here.
-
-  double coefficient (void) const { return coefficient10 (); }
-  int exponent (void) const { return exponent10 (); }
-
-  double coefficient10 (void) const { return c10; }
-  int exponent10 (void) const { return e10; }
-
-  double coefficient2 (void) const { return c2; }
-  int exponent2 (void) const { return e2; }
-
-  double value (void) const;
-
-  friend std::ostream&  operator << (std::ostream& os, const DET& a);
-
-private:
-
-  // Constructed this way, we assume base 2.
-
-  DET (double c, int e)
-    : c2 (c), c10 (0), e2 (e), e10 (0), base2 (true)
-    {
-      initialize10 ();
-    }
-
-  // Original interface had only this constructor and it was assumed
-  // to be base 10, so we are preserving that interface here.
-
-  DET (const double *d)
-    : c2 (0), c10 (d[0]), e2 (0), e10 (static_cast<int> (d[1])), base2 (false)
-    {
-      initialize2 ();
-    }
-
-  void initialize2 (void);
-  void initialize10 (void);
-
-  double c2;
-  double c10;
-
-  int e2;
-  int e10;
-
-  // TRUE means the original values were provided in base 2.
-  bool base2;
-};
-
-#endif
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
--- a/liboctave/fCMatrix.cc	Tue Nov 18 22:55:13 2008 +0100
+++ b/liboctave/fCMatrix.cc	Wed Nov 19 11:23:07 2008 +0100
@@ -38,7 +38,7 @@
 
 #include "Array-util.h"
 #include "fCMatrix.h"
-#include "fCmplxDET.h"
+#include "DET.h"
 #include "fCmplxSCHUR.h"
 #include "fCmplxSVD.h"
 #include "fCmplxCHOL.h"
@@ -1624,33 +1624,13 @@
 	    } 
 	  else 
 	    {
-	      FloatComplex c = 1.0;
-	      int e = 0;
-
+              retval = FloatComplexDET (1.0);
+              
 	      for (octave_idx_type i = 0; i < nc; i++) 
 		{
-		  if (ipvt(i) != (i+1))
-		    c = -c;
-
-		  c *= atmp(i,i);
-
-		  if (c == static_cast<float> (0.0))
-		    break;
-
-		  while (std::abs(c) < 0.5)
-		    {
-		      c *= 2.0;
-		      e--;
-		    }
-
-		  while (std::abs(c) >= 2.0)
-		    {
-		      c /= 2.0;
-		      e++;
-		    }
-		}
-
-	      retval = FloatComplexDET (c, e);
+                  FloatComplex c = atmp(i,i);
+                  retval *= (ipvt(i) != (i+1)) ? -c : c;
+                }
 	    }
 	}
     }
--- a/liboctave/fCMatrix.h	Tue Nov 18 22:55:13 2008 +0100
+++ b/liboctave/fCMatrix.h	Wed Nov 19 11:23:07 2008 +0100
@@ -31,6 +31,7 @@
 #include "mx-defs.h"
 #include "mx-op-defs.h"
 #include "oct-cmplx.h"
+#include "DET.h"
 
 class
 OCTAVE_API
--- a/liboctave/fCmplxDET.cc	Tue Nov 18 22:55:13 2008 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,86 +0,0 @@
-/*
-
-Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2006,
-              2007 John W. Eaton
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <cassert>
-#include <cfloat>
-
-#include "fCmplxDET.h"
-#include "lo-mappers.h"
-#include "lo-math.h"
-#include "oct-cmplx.h"
-
-bool
-FloatComplexDET::value_will_overflow (void) const
-{
-  return base2
-    ? (e2 + 1 > xlog2 (FLT_MAX) ? 1 : 0)
-    : (e10 + 1 > log10 (FLT_MAX) ? 1 : 0);
-}
-
-bool
-FloatComplexDET::value_will_underflow (void) const
-{
-  return base2
-    ? (e2 - 1 < xlog2 (FLT_MIN) ? 1 : 0)
-    : (e10 - 1 < log10 (FLT_MIN) ? 1 : 0);
-}
-
-void
-FloatComplexDET::initialize10 (void)
-{
-  if (c2 != static_cast<float> (0.0))
-    {
-      float etmp = e2 / xlog2 (static_cast<float>(10));
-      e10 = static_cast<int> (xround (etmp));
-      etmp -= e10;
-      c10 = c2 * static_cast<float> (pow (10.0f, etmp));
-    }
-}
-
-void
-FloatComplexDET::initialize2 (void)
-{
-  if (c10 != static_cast<float> (0.0))
-    {
-      float etmp = e10 / log10 (2.0);
-      e2 = static_cast<int> (xround (etmp));
-      etmp -= e2;
-      c2 = c10 * xexp2 (etmp);
-    }
-}
-
-FloatComplex
-FloatComplexDET::value (void) const
-{
-  return base2 ? c2 * xexp2 (static_cast<float>(e2)) : c10 * static_cast<float> (pow (10.0, e10));
-}
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
--- a/liboctave/fCmplxDET.h	Tue Nov 18 22:55:13 2008 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,120 +0,0 @@
-/*
-
-Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
-              2007 John W. Eaton
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#if !defined (octave_FloatComplexDET_h)
-#define octave_FloatComplexDET_h 1
-
-#include <iostream>
-
-#include "oct-cmplx.h"
-
-// FIXME -- we could use templates here; compare with dbleDET.h
-
-class
-OCTAVE_API
-FloatComplexDET
-{
-friend class FloatComplexMatrix;
-
-public:
-
-  FloatComplexDET (void) : c2 (0), c10 (0), e2 (0), e10 (0), base2 (false) { }
-
-  FloatComplexDET (const FloatComplexDET& a)
-    : c2 (a.c2), c10 (a.c10), e2 (a.e2), e10 (a.e10), base2 (a.base2)
-    { }
-
-  FloatComplexDET& operator = (const FloatComplexDET& a)
-    {
-      if (this != &a)
-	{
-	  c2 = a.c2;
-	  e2 = a.e2;
-
-	  c10 = a.c10;
-	  e10 = a.e10;
-
-	  base2 = a.base2;
-	}
-      return *this;
-    }
-
-  bool value_will_overflow (void) const;
-  bool value_will_underflow (void) const;
-
-  // These two functions were originally defined in base 10, so we are
-  // preserving that interface here.
-
-  FloatComplex coefficient (void) const { return coefficient10 (); }
-  int exponent (void) const { return exponent10 (); }
-
-  FloatComplex coefficient10 (void) const { return c10; }
-  int exponent10 (void) const { return e10; }
-
-  FloatComplex coefficient2 (void) const { return c2; }
-  int exponent2 (void) const { return e2; }
-
-  FloatComplex value (void) const;
-
-  friend std::ostream&  operator << (std::ostream& os, const FloatComplexDET& a);
-
-private:
-
-  // Constructed this way, we assume base 2.
-
-  FloatComplexDET (const FloatComplex& c, int e)
-    : c2 (c), c10 (0), e2 (e), e10 (0), base2 (true)
-    {
-      initialize10 ();
-    }
-
-  // Original interface had only this constructor and it was assumed
-  // to be base 10, so we are preserving that interface here.
-
-  FloatComplexDET (const FloatComplex *d)
-    : c2 (0), c10 (d[0]), e2 (0), e10 (static_cast<int> (d[1].real ())),
-      base2 (false)
-    {
-      initialize2 ();
-    }
-
-  void initialize2 (void);
-  void initialize10 (void);
-
-  FloatComplex c2;
-  FloatComplex c10;
-
-  int e2;
-  int e10;
-
-  // TRUE means the original values were provided in base 2.
-  bool base2;
-};
-
-#endif
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
--- a/liboctave/fMatrix.cc	Tue Nov 18 22:55:13 2008 +0100
+++ b/liboctave/fMatrix.cc	Wed Nov 19 11:23:07 2008 +0100
@@ -35,7 +35,7 @@
 #include "Array-util.h"
 #include "byte-swap.h"
 #include "fMatrix.h"
-#include "floatDET.h"
+#include "DET.h"
 #include "floatSCHUR.h"
 #include "floatSVD.h"
 #include "floatCHOL.h"
@@ -1296,33 +1296,13 @@
 	    } 
 	  else 
 	    {
-	      float c = 1.0;
-	      int e = 0;
-
+              retval = FloatDET (1.0);
+              
 	      for (octave_idx_type i = 0; i < nc; i++) 
 		{
-		  if (ipvt(i) != (i+1))
-		    c = -c;
-
-		  c *= atmp(i,i);
-
-		  if (c == 0.0)
-		    break;
-
-		  while (fabs (c) < 0.5)
-		    {
-		      c *= 2.0;
-		      e--;
-		    }
-
-		  while (fabs (c) >= 2.0)
-		    {
-		      c /= 2.0;
-		      e++;
-		    }
-		}
-
-	      retval = FloatDET (c, e);
+                  float c = atmp(i,i);
+                  retval *= (ipvt(i) != (i+1)) ? -c : c;
+                }
 	    }
 	}
     }
--- a/liboctave/fMatrix.h	Tue Nov 18 22:55:13 2008 +0100
+++ b/liboctave/fMatrix.h	Wed Nov 19 11:23:07 2008 +0100
@@ -30,6 +30,7 @@
 
 #include "mx-defs.h"
 #include "mx-op-defs.h"
+#include "DET.h"
 
 class
 OCTAVE_API
--- a/liboctave/floatDET.cc	Tue Nov 18 22:55:13 2008 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,84 +0,0 @@
-/*
-
-Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2006,
-              2007 John W. Eaton
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <cfloat>
-
-#include "floatDET.h"
-#include "lo-mappers.h"
-#include "lo-math.h"
-
-bool
-FloatDET::value_will_overflow (void) const
-{
-  return base2
-    ? (e2 + 1 > xlog2 (FLT_MAX) ? 1 : 0)
-    : (e10 + 1 > log10 (FLT_MAX) ? 1 : 0);
-}
-
-bool
-FloatDET::value_will_underflow (void) const
-{
-  return base2
-    ? (e2 - 1 < xlog2 (FLT_MIN) ? 1 : 0)
-    : (e10 - 1 < log10 (FLT_MIN) ? 1 : 0);
-}
-
-void
-FloatDET::initialize10 (void)
-{
-  if (c2 != 0.0)
-    {
-      float etmp = e2 / xlog2 (static_cast<float>(10));
-      e10 = static_cast<int> (xround (etmp));
-      etmp -= e10;
-      c10 = c2 * powf (10.0, etmp);
-    }
-}
-
-void
-FloatDET::initialize2 (void)
-{
-  if (c10 != 0.0)
-    {
-      float etmp = e10 / log10 (2.0);
-      e2 = static_cast<int> (xround (etmp));
-      etmp -= e2;
-      c2 = c10 * xexp2 (etmp);
-    }
-}
-
-float
-FloatDET::value (void) const
-{
-  return base2 ? c2 * xexp2 (static_cast<float>(e2)) : c10 * powf (10.0, e10);
-}
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
--- a/liboctave/floatDET.h	Tue Nov 18 22:55:13 2008 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,117 +0,0 @@
-/*
-
-Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
-              2007 John W. Eaton
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#if !defined (octave_FloatDET_h)
-#define octave_FloatDET_h 1
-
-#include <iostream>
-
-// FIXME -- we could use templates here; compare with CmplxFloatDET.h
-
-class
-OCTAVE_API
-FloatDET
-{
-friend class FloatMatrix;
-
-public:
-
-  FloatDET (void) : c2 (0), c10 (0), e2 (0), e10 (0), base2 (false) { }
-
-  FloatDET (const FloatDET& a)
-    : c2 (a.c2), c10 (a.c10), e2 (a.e2), e10 (a.e10), base2 (a.base2)
-    { }
-
-  FloatDET& operator = (const FloatDET& a)
-    {
-      if (this != &a)
-	{
-	  c2 = a.c2;
-	  e2 = a.e2;
-
-	  c10 = a.c10;
-	  e10 = a.e10;
-
-	  base2 = a.base2;
-	}
-      return *this;
-    }
-
-  bool value_will_overflow (void) const;
-  bool value_will_underflow (void) const;
-
-  // These two functions were originally defined in base 10, so we are
-  // preserving that interface here.
-
-  float coefficient (void) const { return coefficient10 (); }
-  int exponent (void) const { return exponent10 (); }
-
-  float coefficient10 (void) const { return c10; }
-  int exponent10 (void) const { return e10; }
-
-  float coefficient2 (void) const { return c2; }
-  int exponent2 (void) const { return e2; }
-
-  float value (void) const;
-
-  friend std::ostream&  operator << (std::ostream& os, const FloatDET& a);
-
-private:
-
-  // Constructed this way, we assume base 2.
-
-  FloatDET (float c, int e)
-    : c2 (c), c10 (0), e2 (e), e10 (0), base2 (true)
-    {
-      initialize10 ();
-    }
-
-  // Original interface had only this constructor and it was assumed
-  // to be base 10, so we are preserving that interface here.
-
-  FloatDET (const float *d)
-    : c2 (0), c10 (d[0]), e2 (0), e10 (static_cast<int> (d[1])), base2 (false)
-    {
-      initialize2 ();
-    }
-
-  void initialize2 (void);
-  void initialize10 (void);
-
-  float c2;
-  float c10;
-
-  int e2;
-  int e10;
-
-  // TRUE means the original values were provided in base 2.
-  bool base2;
-};
-
-#endif
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
--- a/liboctave/mx-defs.h	Tue Nov 18 22:55:13 2008 +0100
+++ b/liboctave/mx-defs.h	Wed Nov 19 11:23:07 2008 +0100
@@ -70,11 +70,6 @@
 class FloatCHOL;
 class FloatComplexCHOL;
 
-class DET;
-class ComplexDET;
-class FloatDET;
-class FloatComplexDET;
-
 class EIG;
 
 class HESS;
--- a/liboctave/mx-ext.h	Tue Nov 18 22:55:13 2008 +0100
+++ b/liboctave/mx-ext.h	Wed Nov 19 11:23:07 2008 +0100
@@ -31,10 +31,7 @@
 
 // Result of a Determinant calculation.
 
-#include "dbleDET.h"
-#include "CmplxDET.h"
-#include "floatDET.h"
-#include "fCmplxDET.h"
+#include "DET.h"
 
 // Result of a Cholesky Factorization
 
--- a/src/ChangeLog	Tue Nov 18 22:55:13 2008 +0100
+++ b/src/ChangeLog	Wed Nov 19 11:23:07 2008 +0100
@@ -1,3 +1,7 @@
+2008-11-19  Jaroslav Hajek  <highegg@gmail.com>
+
+	* DLD_FUNCTIONS/det.cc: Include only DET.h.
+
 2008-11-17  John W. Eaton  <jwe@octave.org>
 
 	* load-path.cc (load_path::dir_info::update): Simplify previous
--- a/src/DLD-FUNCTIONS/det.cc	Tue Nov 18 22:55:13 2008 +0100
+++ b/src/DLD-FUNCTIONS/det.cc	Wed Nov 19 11:23:07 2008 +0100
@@ -25,10 +25,7 @@
 #include <config.h>
 #endif
 
-#include "CmplxDET.h"
-#include "dbleDET.h"
-#include "fCmplxDET.h"
-#include "floatDET.h"
+#include "DET.h"
 
 #include "defun-dld.h"
 #include "error.h"