# HG changeset patch
# User jwe
# Date 1141270801 0
# Node ID 4b45b2bcda893511f2101a4f09203eb7156c37f2
# Parent  92f8b2723c8cfe86b339a2f117bd953cea6f72dc
[project @ 2006-03-02 03:40:00 by jwe]

diff --git a/liboctave/CMatrix.cc b/liboctave/CMatrix.cc
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -1439,10 +1439,7 @@
 
   if (nr == 0 || nc == 0)
     {
-      Complex d[2];
-      d[0] = 1.0;
-      d[1] = 0.0;
-      retval = ComplexDET (d);
+      retval = ComplexDET (1.0, 0);
     }
   else
     {
@@ -1500,24 +1497,33 @@
 		} 
 	      else 
 		{
-		  Complex d[2] = { 1., 0.};
-		  for (octave_idx_type i=0; i<nc; i++) 
+		  Complex c = 1.0;
+		  int e = 0;
+
+		  for (octave_idx_type i = 0; i < nc; i++) 
 		    {
-		      if (ipvt(i) != (i+1)) d[0] = -d[0];
-		      d[0] = d[0] * atmp(i,i);
-		      if (d[0] == 0.) break;
-		      while (std::abs(d[0]) < 1.) 
+		      if (ipvt(i) != (i+1))
+			c = -c;
+
+		      c *= atmp(i,i);
+
+		      if (c == 0.0)
+			break;
+
+		      while (std::abs(c) < 0.5)
 			{
-			  d[0] = 10. * d[0];
-			  d[1] = d[1] - 1.0;
+			  c *= 2.0;
+			  e--;
 			}
-		      while (std::abs(d[0]) >= 10.) 
+
+		      while (std::abs(c) >= 2.0)
 			{
-			  d[0] = 0.1 * d[0];
-			  d[1] = d[1] + 1.0;
+			  c /= 2.0;
+			  e++;
 			}
 		    }
-		  retval = ComplexDET (d);
+
+		  retval = ComplexDET (c, e);
 		}
 	    }
 	}
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,34 @@
+2006-03-01  John W. Eaton  <jwe@octave.org>
+
+	* CMatrix.cc (ComplexMatrix::determinant):
+	Scale result by factors of 2, not 10.
+	* dMatrix.cc (Matrix::determinant): Likewise.
+
+	* dbleDET.h (DET::DET): Use initializer list.
+	(DET::coefficient2, DET::coefficient10, DET::exponent2,
+	DET::exponent10): New functions.
+	(DET::det): Delete.
+	(DET::c2, DET::c10, DET::e2, DET::e10, DET::base2): New data members.
+	Store value internally with double and int instead of 2-element
+	double vector.
+	(DET::initialize2, DET::initialize10): Provide decls.
+	* dbleDET.cc (DET::value_will_overflow,	DET::value_will_underflow): 
+	Return bool value, not int.
+	(DET::initialize2, DET::initialize10): New functions.
+
+	* CmplxDET.h (ComplexDET::ComplexDET): Use initializer list.
+	(ComplexDET::coefficient2, ComplexDET::coefficient10,
+	ComplexDET::exponent2, ComplexDET::exponent10): New functions.
+	(ComplexDET::det): Delete.
+	(ComplexDET::c2, ComplexDET::c10, ComplexDET::e2, ComplexDET::e10,
+	ComplexDET::base2): New data members.
+	Store value internally with Complex and int instead of 2-element
+	Complex vector.
+	(ComplexDET::initialize2, ComplexDET::initialize10): Provide decls.
+	* dbleComplexDET.cc (ComplexDET::value_will_overflow,
+	ComplexDET::value_will_underflow): Return bool value, not int.
+	(ComplexDET::initialize2, ComplexDET::initialize10): New functions.
+
 2006-02-24  John W. Eaton  <jwe@octave.org>
 
 	* Array.cc (assignN): Clear index before reshaping.
diff --git a/liboctave/CmplxDET.cc b/liboctave/CmplxDET.cc
--- a/liboctave/CmplxDET.cc
+++ b/liboctave/CmplxDET.cc
@@ -25,40 +25,58 @@
 #include <config.h>
 #endif
 
+#include <cassert>
 #include <cfloat>
 #include <cmath>
 
 #include "CmplxDET.h"
+#include "lo-mappers.h"
 #include "oct-cmplx.h"
 
-int
+bool
 ComplexDET::value_will_overflow (void) const
 {
-  return det[1].real () + 1 > log10 (DBL_MAX) ? 1 : 0;
+  return base2
+    ? (e2 + 1 > log2 (DBL_MAX) ? 1 : 0)
+    : (e10 + 1 > log10 (DBL_MAX) ? 1 : 0);
 }
 
-int
+bool
 ComplexDET::value_will_underflow (void) const
 {
-  return det[1].real () - 1 < log10 (DBL_MIN) ? 1 : 0;
+  return base2
+    ? (e2 - 1 < log2 (DBL_MIN) ? 1 : 0)
+    : (e10 - 1 < log10 (DBL_MIN) ? 1 : 0);
 }
 
-Complex
-ComplexDET::coefficient (void) const
+void
+ComplexDET::initialize10 (void)
 {
-  return det[0];
+  if (c2 != 0.0)
+    {
+      double etmp = e2 / log2 (10);
+      e10 = static_cast<int> (round (etmp));
+      etmp -= e10;
+      c10 = c2 * pow (10.0, etmp);
+    }
 }
 
-int
-ComplexDET::exponent (void) const
+void
+ComplexDET::initialize2 (void)
 {
-  return (int) (det[1].real ());
+  if (c10 != 0.0)
+    {
+      double etmp = e10 / log10 (2);
+      e2 = static_cast<int> (round (etmp));
+      etmp -= e2;
+      c2 = c10 * exp2 (etmp);
+    }
 }
 
 Complex
 ComplexDET::value (void) const
 {
-  return det[0] * pow (10.0, det[1].real ());
+  return base2 ? c2 * exp2 (e2) : c10 * pow (10.0, e10);
 }
 
 /*
diff --git a/liboctave/CmplxDET.h b/liboctave/CmplxDET.h
--- a/liboctave/CmplxDET.h
+++ b/liboctave/CmplxDET.h
@@ -28,6 +28,8 @@
 
 #include "oct-cmplx.h"
 
+// XXX FIXME XXX -- we could use templates here; compare with dbleDET.h
+
 class
 ComplexDET
 {
@@ -36,30 +38,41 @@
 
 public:
 
-  ComplexDET (void) { }
+  ComplexDET (void) : c2 (0), c10 (0), e2 (0), e10 (0), base2 (false) { }
 
   ComplexDET (const ComplexDET& a)
-    {
-      det[0] = a.det[0];
-      det[1] = a.det[1];
-    }
+    : c2 (a.c2), c10 (a.c10), e2 (a.e2), e10 (a.e10), base2 (a.base2)
+    { }
 
   ComplexDET& operator = (const ComplexDET& a)
     {
       if (this != &a)
 	{
-	  det[0] = a.det[0];
-	  det[1] = a.det[1];
+	  c2 = a.c2;
+	  e2 = a.e2;
+
+	  c10 = a.c10;
+	  e10 = a.e10;
+
+	  base2 = a.base2;
 	}
       return *this;
     }
 
-  int value_will_overflow (void) const;
-  int value_will_underflow (void) const;
+  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;
+  Complex coefficient (void) const { return coefficient10 (); }
+  int exponent (void) const { return exponent10 (); }
 
-  int exponent (void) const;
+  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;
 
@@ -67,13 +80,35 @@
 
 private:
 
-  ComplexDET (const Complex *d)
+  // Constructed this way, we assume base 2.
+
+  ComplexDET (const Complex& c, int e)
+    : c2 (c), c10 (0), e2 (e), e10 (0), base2 (true)
     {
-      det[0] = d[0];
-      det[1] = d[1];
+      initialize10 ();
     }
 
-  Complex det [2];
+  // 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
diff --git a/liboctave/dMatrix.cc b/liboctave/dMatrix.cc
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -1102,10 +1102,7 @@
 
   if (nr == 0 || nc == 0)
     {
-      double d[2];
-      d[0] = 1.0;
-      d[1] = 0.0;
-      retval = DET (d);
+      retval = DET (1.0, 0);
     }
   else
     {
@@ -1163,24 +1160,33 @@
 		} 
 	      else 
 		{
-		  double d[2] = { 1., 0.};
-		  for (octave_idx_type i=0; i<nc; i++) 
+		  double c = 1.0;
+		  int e = 0;
+
+		  for (octave_idx_type i = 0; i < nc; i++) 
 		    {
-		      if (ipvt(i) != (i+1)) d[0] = -d[0];
-		      d[0] *= atmp(i,i);
-		      if (d[0] == 0.) break;
-		      while (fabs(d[0]) < 1.) 
+		      if (ipvt(i) != (i+1))
+			c = -c;
+
+		      c *= atmp(i,i);
+
+		      if (c == 0.0)
+			break;
+
+		      while (fabs (c) < 0.5)
 			{
-			  d[0] = 10. * d[0];
-			  d[1] = d[1] - 1.0;
+			  c *= 2.0;
+			  e--;
 			}
-		      while (fabs(d[0]) >= 10.) 
+
+		      while (fabs (c) >= 2.0)
 			{
-			  d[0] = 0.1 * d[0];
-			  d[1] = d[1] + 1.0;
+			  c /= 2.0;
+			  e++;
 			}
 		    }
-		  retval = DET (d);
+
+		  retval = DET (c, e);
 		}
 	    }
 	}
diff --git a/liboctave/dbleDET.cc b/liboctave/dbleDET.cc
--- a/liboctave/dbleDET.cc
+++ b/liboctave/dbleDET.cc
@@ -30,34 +30,50 @@
 
 #include "dbleDET.h"
 
-int
+bool
 DET::value_will_overflow (void) const
 {
-  return det[1] + 1 > log10 (DBL_MAX) ? 1 : 0;
+  return base2
+    ? (e2 + 1 > log2 (DBL_MAX) ? 1 : 0)
+    : (e10 + 1 > log10 (DBL_MAX) ? 1 : 0);
 }
 
-int
+bool
 DET::value_will_underflow (void) const
 {
-  return det[1] - 1 < log10 (DBL_MIN) ? 1 : 0;
+  return base2
+    ? (e2 - 1 < log2 (DBL_MIN) ? 1 : 0)
+    : (e10 - 1 < log10 (DBL_MIN) ? 1 : 0);
 }
 
-double
-DET::coefficient (void) const
+void
+DET::initialize10 (void)
 {
-  return det[0];
+  if (c2 != 0.0)
+    {
+      double etmp = e2 / log2 (10);
+      e10 = static_cast<int> (round (etmp));
+      etmp -= e10;
+      c10 = c2 * pow (10.0, etmp);
+    }
 }
 
-int
-DET::exponent (void) const
+void
+DET::initialize2 (void)
 {
-  return (int) det[1];
+  if (c10 != 0.0)
+    {
+      double etmp = e10 / log10 (2);
+      e2 = static_cast<int> (round (etmp));
+      etmp -= e2;
+      c2 = c10 * exp2 (etmp);
+    }
 }
 
 double
 DET::value (void) const
 {
-  return det[0] * pow (10.0, det[1]);
+  return base2 ? c2 * exp2 (e2) : c10 * pow (10.0, e10);
 }
 
 /*
diff --git a/liboctave/dbleDET.h b/liboctave/dbleDET.h
--- a/liboctave/dbleDET.h
+++ b/liboctave/dbleDET.h
@@ -26,6 +26,8 @@
 
 #include <iostream>
 
+// XXX FIXME XXX -- we could use templates here; compare with CmplxDET.h
+
 class
 DET
 {
@@ -34,30 +36,41 @@
 
 public:
 
-  DET (void) { }
+  DET (void) : c2 (0), c10 (0), e2 (0), e10 (0), base2 (false) { }
 
   DET (const DET& a)
-    {
-      det[0] = a.det[0];
-      det[1] = a.det[1];
-    }
+    : c2 (a.c2), c10 (a.c10), e2 (a.e2), e10 (a.e10), base2 (a.base2)
+    { }
 
   DET& operator = (const DET& a)
     {
       if (this != &a)
 	{
-	  det[0] = a.det[0];
-	  det[1] = a.det[1];
+	  c2 = a.c2;
+	  e2 = a.e2;
+
+	  c10 = a.c10;
+	  e10 = a.e10;
+
+	  base2 = a.base2;
 	}
       return *this;
     }
 
-  int value_will_overflow (void) const;
-  int value_will_underflow (void) const;
+  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;
+  double coefficient (void) const { return coefficient10 (); }
+  int exponent (void) const { return exponent10 (); }
 
-  int exponent (void) const;
+  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;
 
@@ -65,13 +78,34 @@
 
 private:
 
-  DET (const double *d)
+  // Constructed this way, we assume base 2.
+
+  DET (double c, int e)
+    : c2 (c), c10 (0), e2 (e), e10 (0), base2 (true)
     {
-      det[0] = d[0];
-      det[1] = d[1];
+      initialize10 ();
     }
 
-  double det [2];
+  // 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