[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Possible bug in det (was: Possible loss of accuracy)
From: |
Stephen Montgomery-Smith |
Subject: |
Re: Possible bug in det (was: Possible loss of accuracy) |
Date: |
Thu, 16 May 2013 12:34:17 -0500 |
User-agent: |
Mozilla/5.0 (X11; Linux i686; rv:17.0) Gecko/20130510 Thunderbird/17.0.6 |
On 05/16/2013 12:26 PM, Ed Meyer wrote:
>
>
> On Thu, May 16, 2013 at 8:30 AM, Stephen Montgomery-Smith
> <address@hidden <mailto:address@hidden>> wrote:
>
> On 05/16/2013 08:47 AM, Jordi GutiƩrrez Hermoso wrote:
>
> > The actual determinant (base_det::value) is then computed as the
> > result of ldexp () after taking the significand and exponent from
> > frexp (). But the problem is that here the significand gets very small
> > and the exponent gets very big, so ldexp () ends up reporting
> > infinity.
> >
> > Stepping through the code, this does indeed look like a bug. I don't
> > understand why it's doing this at all instead of just ordinary
> > floating point multiplication. I'm guessing it's trying to work around
> > problems when losing precision due to multiplying numbers that are of
> > very different magnitude. Please open a bug report, and perhaps
> > someone else can figure out what the problem is.
>
> Maybe they are trying to get around the problem where the product of the
> diagonal entries is a fairly reasonably sized number, but the diagonal
> entries consists of very large and very small numbers, so that if you
> multiply the diagonal entries in a certain order, you will get an
> underflow or an overflow.
>
> Straight multiplying is not going to fix this problem.
>
> I think a better approach is to separate the mantissa from the exponent,
> and then keep a running total of the products of the mantissa and the
> sum of the exponents, BUT also if the product of the mantissa gets too
> large or too small where there is a danger of overflow or underflow, the
> mantissa should be appropriately modified.
>
> So change the code in:
>
> http://hg.savannah.gnu.org/hgweb/octave/file/7f6f0b3f7369/liboctave/numeric/DET.h#l70
>
> to:
>
> c2 *= xlog2 (t, e);
> e2 += e;
> c2 = xlog2 (c2, e);
> e2 += e;
>
> or:
>
> c2 = xlog2 (c2*t, e);
> e2 += e;
>
> I think the first option might be marginally safer if t is very close to
> DOUBLE_MAX or DOUBLE_MIN.
>
>
> _______________________________________________
> Help-octave mailing list
> address@hidden <mailto:address@hidden>
> https://mailman.cae.wisc.edu/listinfo/help-octave
>
>
> More efficient and probably more accurate is the way the old linpack
> did it, by keeping two numbers, one lying between 1 and 10, the other
> the exponent of ten:
>
> http://www.netlib.org/linpack/sgedi.f
>
> then you could flag a retval that will overflow
But the algorithm I was suggesting is the same, except 10 is replaced by
2. And powers of 2 are so much easier to handle than powers of 10.