Re: Possible bug in det (was: Possible loss of accuracy)

[Stephen Montgomery-Smith]

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.