[bug #26190] dec2base test fails due to floating point error

[Thorsten Meyer]

URL:
  <http://savannah.gnu.org/bugs/?26190>

                 Summary: dec2base test fails due to floating point error
                 Project: GNU Octave
            Submitted by: tmeyier
            Submitted on: Wed 15 Apr 2009 05:16:36 PM UTC
                Category: Libraries
                Severity: 3 - Normal
              Item Group: Incorrect Result
                  Status: Need Info
             Assigned to: tmeyier
         Originator Name: 
        Originator Email: 
             Open/Closed: Open
         Discussion Lock: Any
                 Release: dev
        Operating System: GNU/Linux

    _______________________________________________________

Details:

with the present tip (bb62bc406ea7) (and also with some other recent build,
if I remember correctly) I get two fails for the tests in dec2base.
This simple example demonstrates the problem:
  octave:160> dec2base(999, 10)
  ans = 998

Looking into the code, the following lines seem to be responsible for the
bug:
  power = ones (length (n), 1) * (base .^ (max_len-1 : -1 : 0));
  n = n(:) * ones (1, max_len);
  digits = floor (rem (n, base*power) ./ power);

power is assigned a vector of powers of the base. However, due to floating
point round off error, the numbers in power are not all integers. For the
above example:
octave:161> max_len=3;
octave:162> base=10;
octave:163> n=999;
octave:167> max_len = round (log (max (max (n), 1)) ./ log (base)) + 1
octave:168> power = ones (length (n), 1) * (base .^ (max_len-1 : -1 : 0))
power =

   1000.0000    100.0000     10.0000      1.0000

octave:170> power-1
ans =

   9.9900e+02   9.9000e+01   9.0000e+00   2.2204e-16

octave:171> n = n(:) * ones (1, max_len);
octave:172>  digits = floor (rem (n, base*power) ./ power)
digits =

   0   9   9   8

power(4) is a bit larger than 1, hence the remainder is a bit too small
leading to the result being rounded down to 8.

Changing the definition of power to
 power = ones (length (n), 1) * round (base .^ (max_len-1 : -1 : 0));

seems to fix the problem. 

Still, before making a changeset I would like to ask the floating point
experts among you if that is the right thing to do.

regards

Thorsten






    _______________________________________________________

Reply to this item at:

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

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