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

[Clemens Buchacher]

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;
     }