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

[louis scott]

Clemens,

I meant to submit a bug report, as the det problem was originally sent by me,  I apologize as I have yet to do so.

See Jordi's and my comments, given that det does an LU decomposition first, then all you need is det(U)

and use http://hg.savannah.gnu.org/hgweb/octave/file/7f6f0b3f7369/liboctave/array/dMatrix.cc#l1263

This is I think what I saw as the answer since no matter what rounding comes to the off diagonals, the zeros in the lower diagonals will kill off the unwanted terms.
N = 2100;
A = rand(N);
det(A*inv(A)) % The correct answer is the det(eye(N)) = 1

I did not focus on the overloaded *= , this may cure some problems in addition, but in the example I shared the above works and only requires an additional N multiplications to compute rather than the full matrix det computations. 

Thank you for so quickly focusing on this. 

Best,

Louis

On Sun, May 19, 2013 at 1:59 PM, Clemens Buchacher <address@hidden> wrote:

Hi,

After testing I am planning on submitting this patch. Is the description
ok with you? Any other comments?

Clemens
---

# HG changeset patch
# User Clemens Buchacher <address@hidden>
# Date 1368966676 -7200
#      Sun May 19 14:31:16 2013 +0200
# Node ID 784e1a2ae039aa0fad8f5d8d9fe0c9e6e3821012
# Parent  a463f6b18f7de3d79d511ddb4324b0c6a1074410
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
Reported-by: Marco Caliari <address@hidden>
Fix-proposed-by: Stephen Montgomery-Smith <address@hidden>

diff -r a463f6b18f7d -r 784e1a2ae039 liboctave/numeric/DET.h
--- a/liboctave/numeric/DET.h   Sat May 11 11:58:47 2013 +0200
+++ b/liboctave/numeric/DET.h   Sun May 19 14:31:16 2013 +0200
@@ -71,7 +71,8 @@
   void operator *= (T t)
     {
       int e;
-      c2 *= xlog2 (t, e);
+      c2 *= t;
+      c2 = xlog2(c2, e);
       e2 += e;

     }

_______________________________________________
Help-octave mailing list
address@hidden
https://mailman.cae.wisc.edu/listinfo/help-octave