Mercurial > hg > octave-nkf
changeset 5634:4b45b2bcda89
[project @ 2006-03-02 03:40:00 by jwe]
author | jwe |
---|---|
date | Thu, 02 Mar 2006 03:40:01 +0000 |
parents | 92f8b2723c8c |
children | fce700d5bb5f |
files | liboctave/CMatrix.cc liboctave/ChangeLog liboctave/CmplxDET.cc liboctave/CmplxDET.h liboctave/dMatrix.cc liboctave/dbleDET.cc liboctave/dbleDET.h |
diffstat | 7 files changed, 230 insertions(+), 84 deletions(-) [+] |
line wrap: on
line diff
--- 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); } } }
--- 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.
--- 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); } /*
--- 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
--- 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); } } }
--- 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); } /*
--- 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