[Octave-bug-tracker] [bug #39014] Wrong determinant for some (large) matrices

[Clemens Buchacher]

Follow-up Comment #1, bug #39014 (project octave):

For the determinant we make some effort to deal with large intermediate
results. We normalize factors using frexp and keep track of the exponent
separately, such that we are not limited by the double precision exponent (+/-
1023). Let's assume all factors are close to 1. For those smaller than 1,
frexp returns an exponent of 0 and leaves the mantissa untouched. For those
equal to or larger than 1, frexp returns an exponent of 1 and divides the
mantissa by 2. So for those latter factors we have a mantissa close to 0.5,
and we keep multiplying with 0.5 until we do hit the double precision limit in
the mantissa:

> det(matrix_type(full(diag(ones(1, 2000)-1e-10)), "full"))
ans =  1.00000
> det(matrix_type(full(diag(ones(1, 2000))), "full"))
ans = NaN

Breakpoint 6, Matrix::determinant (address@hidden, 
    mattype=..., address@hidden: 0, address@hidden: 1, 
    address@hidden) at ../../liboctave/array/dMatrix.cc:1362
1362                          retval *= (ipvt(i) != (i+1)) ? -c : c;
3: i = 1066
2: retval = {c2 = 6.3240402667679558e-322, e2 = 1067}
(gdb) 
Will ignore next 9 crossings of breakpoint 6.  Continuing.

Breakpoint 6, Matrix::determinant (address@hidden, 
    mattype=..., address@hidden: 0, address@hidden: 1, 
    address@hidden) at ../../liboctave/array/dMatrix.cc:1362
1362                          retval *= (ipvt(i) != (i+1)) ? -c : c;
3: i = 1076
2: retval = {c2 = 0, e2 = 1077}

I wonder why the limit on my machine appears to be 2^-1074 ~ 1e-325.

So considering the effort we make already, it seems that we should be doing
better here. How about renormalizing retval every 100 iterations or so?

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?39014>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/