Mon 20 May 2013 09:07:18 PM UTC, comment #5:
Clemens Buchacher wrote:
> on yet another octave mailing list.
I am not sure if you're suggesting that we have too many Octave mailing lists, but if so, you should rejoice about the situation getting better. We had three, help, dev, and OctaveForge, but we recently folded OctaveForge into the dev list.
There are now only two Octave lists: help and dev.

Mon 20 May 2013 01:48:52 PM UTC, comment #4:
Attaching det.patch with the following description:
normalize multiplication result instead of operand
The operand of a multiplication with base_det is normalized to c * 2^e,
such that 0.5 <= c < 1. The mantissa is multiplied by c, and the
exponent is updated separately. After each multiplication, the mantissa
is therefore smaller by a up to a factor 2. Despite the normalization of
the operands, the mantissa can therefore become very small. Eventually,
each multiplication suffers from severe loss of precision and the
mantissa can even get rounded to 0.
Instead, normalize the mantissa after multiplication with the operand.
This makes the operation insensitive to the number of operands.
The change may degrade accuracy in case of operands close to the maximum
or minimum representable floating point numbers, but that is not a
problem base_det tries to solve.
Bug report: https://savannah.gnu.org/bugs/?39014
Reportedby: Marco Caliari <marco.caliari@univr.it>
Fixproposedby: Stephen MontgomerySmith <stephen@missouri.edu>
(file #28118)

Sat 18 May 2013 11:03:09 AM UTC, comment #2:
Err, I veered off course a bit here. But I suspect the original issue is quite similar:
octave3.7.3+:52> det(matrix_type(full(diag(ones(1,2000)+1e10)), "full"))
ans = Inf
In this case, we do not end up with c2=0, but we end up with a very small mantissa, which due to loss of precision does not grow as quickly as the exponent. So the final result ends up as retval = c2 * 2^e2 = 2^1074 * 2^2000, which I suppose gets converted to Inf somewhere (maybe we could do better there too, and get 2^928, but that would not save our day).
Breakpoint 6, Matrix::determinant (this=this@entry=0x7fffffffd2e0,
mattype=..., info=@0x7fffffffd250: 0, rcon=@0x7fffffffd280: 1,
calc_cond=calc_cond@entry=1) at ../../liboctave/array/dMatrix.cc:1362
1362 retval *= (ipvt(i) != (i+1)) ? c : c;
4: c = 1.0000000001
3: i = 1821
2: retval = {c2 = 4.9406564584124654e324, e2 = 1822}

Sat 18 May 2013 10:49:51 AM UTC, comment #1:
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)1e10)), "full"))
ans = 1.00000
> det(matrix_type(full(diag(ones(1, 2000))), "full"))
ans = NaN
Breakpoint 6, Matrix::determinant (this=this@entry=0x7fffffffd2e0,
mattype=..., info=@0x7fffffffd250: 0, rcon=@0x7fffffffd280: 1,
calc_cond=calc_cond@entry=1) at ../../liboctave/array/dMatrix.cc:1362
1362 retval *= (ipvt(i) != (i+1)) ? c : c;
3: i = 1066
2: retval = {c2 = 6.3240402667679558e322, e2 = 1067}
(gdb)
Will ignore next 9 crossings of breakpoint 6. Continuing.
Breakpoint 6, Matrix::determinant (this=this@entry=0x7fffffffd2e0,
mattype=..., info=@0x7fffffffd250: 0, rcon=@0x7fffffffd280: 1,
calc_cond=calc_cond@entry=1) 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 ~ 1e325.
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?

Fri 17 May 2013 01:37:41 PM UTC, original submission:
Dear all,
as reported in the thread https://mailman.cae.wisc.edu/pipermail/helpoctave/2013May/058753.html, the following code
quite often gives Inf instead of 1. The problem seems to be in the way the product of the diagonal elements of the R (or U) factor of the matrix is done. I don't know if prod is sufficiently robust with respect to the product of numbers very different in magnitude. I see
but
Suppose prod is fixed if needed. Isn't then it possible to rewrite det function with a simple m file basically doing
By the way, the actual implementation of prod allows the script above to work fine on A*inv(A).
Cheers,
Marco
