diff options
author | Gael Guennebaud <g.gael@free.fr> | 2009-07-16 14:20:36 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2009-07-16 14:20:36 +0200 |
commit | 525da6a464c897d0fe4e401a65851bcd7567fc5a (patch) | |
tree | 29567b1ed5d365576a094dc3919776d96d1aac1d /Eigen | |
parent | 65fc70b75039a5cdfc5df67f62d38b317196293b (diff) |
bugfix in blueNorm
Diffstat (limited to 'Eigen')
-rw-r--r-- | Eigen/src/Core/Dot.h | 39 | ||||
-rw-r--r-- | Eigen/src/Core/NumTraits.h | 8 |
2 files changed, 10 insertions, 37 deletions
diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h index c5f2e8505..d97f5837f 100644 --- a/Eigen/src/Core/Dot.h +++ b/Eigen/src/Core/Dot.h @@ -304,29 +304,6 @@ MatrixBase<Derived>::stableNorm() const return this->cwise().abs().redux(ei_scalar_hypot_op<RealScalar>()); } -/** \internal Computes ibeta^iexp by binary expansion of iexp, - * exact if ibeta is the machine base */ -template<typename T> inline T bexp(int ibeta, int iexp) -{ - T tbeta = T(ibeta); - T res = 1.0; - int n = iexp; - if (n<0) - { - n = - n; - tbeta = 1.0/tbeta; - } - for(;;) - { - if ((n % 2)==0) - res = res * tbeta; - n = n/2; - if (n==0) return res; - tbeta = tbeta*tbeta; - } - return res; -} - /** \returns the \em l2 norm of \c *this using the Blue's algorithm. * A Portable Fortran Program to Find the Euclidean Norm of a Vector, * ACM TOMS, Vol 4, Issue 1, 1978. @@ -337,7 +314,7 @@ template<typename Derived> inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real MatrixBase<Derived>::blueNorm() const { - static int nmax; + static int nmax = -1; static Scalar b1, b2, s1m, s2m, overfl, rbig, relerr; int n; Scalar ax, abig, amed, asml; @@ -355,8 +332,8 @@ MatrixBase<Derived>::blueNorm() const // are used. For any specific computer, each of the assignment // statements can be replaced nbig = std::numeric_limits<int>::max(); // largest integer - ibeta = NumTraits<Scalar>::Base; // base for floating-point numbers - it = NumTraits<Scalar>::Mantissa; // number of base-beta digits in mantissa + ibeta = std::numeric_limits<Scalar>::radix; //NumTraits<Scalar>::Base; // base for floating-point numbers + it = std::numeric_limits<Scalar>::digits; //NumTraits<Scalar>::Mantissa; // number of base-beta digits in mantissa iemin = std::numeric_limits<Scalar>::min_exponent; // minimum exponent iemax = std::numeric_limits<Scalar>::max_exponent; // maximum exponent rbig = std::numeric_limits<Scalar>::max(); // largest floating-point number @@ -368,17 +345,17 @@ MatrixBase<Derived>::blueNorm() const ei_assert(false && "the algorithm cannot be guaranteed on this computer"); } iexp = -((1-iemin)/2); - b1 = bexp<Scalar>(ibeta, iexp); // lower boundary of midrange + b1 = std::pow(ibeta, iexp); // lower boundary of midrange iexp = (iemax + 1 - it)/2; - b2 = bexp<Scalar>(ibeta,iexp); // upper boundary of midrange + b2 = std::pow(ibeta,iexp); // upper boundary of midrange iexp = (2-iemin)/2; - s1m = bexp<Scalar>(ibeta,iexp); // scaling factor for lower range + s1m = std::pow(ibeta,iexp); // scaling factor for lower range iexp = - ((iemax+it)/2); - s2m = bexp<Scalar>(ibeta,iexp); // scaling factor for upper range + s2m = std::pow(ibeta,iexp); // scaling factor for upper range overfl = rbig*s2m; // overfow boundary for abig - eps = bexp<Scalar>(ibeta, 1-it); + eps = std::pow(ibeta, 1-it); relerr = ei_sqrt(eps); // tolerance for neglecting asml abig = 1.0/eps - 1.0; if (Scalar(nbig)>abig) nmax = abig; // largest safe n diff --git a/Eigen/src/Core/NumTraits.h b/Eigen/src/Core/NumTraits.h index dec07a036..24afe54c5 100644 --- a/Eigen/src/Core/NumTraits.h +++ b/Eigen/src/Core/NumTraits.h @@ -70,9 +70,7 @@ template<> struct NumTraits<float> HasFloatingPoint = 1, ReadCost = 1, AddCost = 1, - MulCost = 1, - Base = 2, - Mantissa = 23 + MulCost = 1 }; }; @@ -85,9 +83,7 @@ template<> struct NumTraits<double> HasFloatingPoint = 1, ReadCost = 1, AddCost = 1, - MulCost = 1, - Base = 2, - Mantissa = 52 + MulCost = 1 }; }; |