On Thu, May 16, 2013 at 8:30 AM, Stephen Montgomery-Smith
<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.