bugGNU Octave - Bugs: bug #31579, Floating point mod function does...

 
 

bug #31579: Floating point mod function does not match Matlab

Submitter:  Rik <rik5>
Submitted:  Sun 07 Nov 2010 06:51:21 PM UTC
   
 
Category:  Interpreter Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  kodo Open/Closed:  * Closed
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Tue 09 Nov 2010 07:08:28 PM UTC, comment #5: 

I checked in the following change for the copysign problem:

http://hg.savannah.gnu.org/hgweb/octave/rev/009d16b010fa

Yes, I think we should report the problem.  The current behavior of fmod does not seem good to me, but maybe there is some reason for it?  Though if there is, I can't really see it.

After doing some debugging, it looks like fmod is ultimately computed by the following function in glibc, at least on my system:

http://sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/ieee754/dbl-64/e_fmod.c;h=2ce613574a0d7aad233dc2ba77036309f45681fd;hb=HEAD

Yowza.  That level of bit twiddling is far beyond my ability to comprehend, much less fix.

In any case, even if fmod is fixed, I think we still need our versions of rem and mod so that we can be compatible with Matlab's "round if within eps" feature.

John W. Eaton <jwe>
Group administrator
Tue 09 Nov 2010 06:20:41 PM UTC, comment #4: 

Do we bother to inform the glibc folks that something is odd with their implementation of fmod?

I used some extended features of C99 to set the rounding mode used by ordinary math operations and it appears that the fmod function is using floor in place of trunc, even though the documentation for fmod says that trunc should be used.


Compile with extra option '-std=c++0x'

#include <iostream>
#include <cmath>
#include <cfenv>

using namespace std;

int
main (void)
{
  double x = 2.1;
  double y = 0.1;

  int rnd_mode = fegetround();
  cout << "Current round mode is " << rnd_mode << endl;
  fesetround(FE_DOWNWARD);
  rnd_mode = fegetround();
  cout << "Current round mode is " << rnd_mode << endl;

  cout << "x = " << x << endl;
  cout << "y = " << y << endl;
  cout << "floor (x / y) = " << floor (x / y) << endl;
  cout << "x - y * floor (x / y) = " << x - y * floor (x / y) << endl;
  cout << "fmod (x, y) = " << fmod (x, y) << endl;
  cout << "(int) (x / y) = " << static_cast<int> (x / y) << endl;
  cout << "x - y * (int) (x / y) = " << x - y * static_cast<int> (x / y) << endl;
  cout << "fmodf (float (x), float (y)) = " << fmodf (float (x), float (y)) << endl;
  cout << "remainder (x, y) = " << remainder (-x, y) << endl;

  return 0;
}


The results are


Current round mode is 0
Current round mode is 1024
x = 2.1
y = 0.1
floor (x / y) = 20
x - y * floor (x / y) = 0.1
fmod (x, y) = 0.1
(int) (x / y) = 20
x - y * (int) (x / y) = 0.1
fmodf (float (x), float (y)) = 2.98023e-08
remainder (x, y) = 2.77556e-17


Rik <rik5>
Group administrator
Tue 09 Nov 2010 06:14:42 PM UTC, comment #3: 

The new mod/rem routines always copy the sign bit.  This does no harm when the remainder is zero, but does look a bit odd.  For example:


mod (2,-1)
ans = -0

whereas in 3.2.4
mod (2,-1)
ans = 0


The change below would fix it, but it isn't critical to implement it.


original:
if (x != y && y != 0)
  retval = copysignf (retval, y);
new:
if (x != y && y != 0 && retval !=0)
  retval = copysignf (retval, y);


Rik <rik5>
Group administrator
Tue 09 Nov 2010 06:17:40 AM UTC, comment #2: 

I checked in the following change which I think makes mod and rem behave approximately the same as in Matlab:

  http://hg.savannah.gnu.org/hgweb/octave/rev/94d9d412a2a0

I say approximately because I can't tell from the description in the Matlab docs precisely how "within roundoff error of an integer" is computed, so I took my best guess.

With this change, mod (2.1, 0.1) returns 0 and I get 2500 for q in the test case from the original report, which I think is bug #31512, not 31514.

As I think my changeset fixes this problem, I'm closing this report.

John W. Eaton <jwe>
Group administrator
Sun 07 Nov 2010 07:58:35 PM UTC, comment #1: 

Octave 3.2.4 gets the expected answer.  It uses the script mod.m which directly calculates:

r = x - y .* floor (x ./ y);


In contrast, the development version of Octave now relies on C libraries for the mod and rem functions.  At least on two systems, JWE's and mine, the fmod function yields strange results.  Attached is C++ code to calculate the mod operation directly, and via the library call to fmod.  The results are different and don't depend on whether the variable is a float or double.

I also re-coded the test into just C, to make sure there were no C++ calling C library function issues, and fmod still produces the wrong result.  This is with glibc 2.10.1 on a 64-bit kernel.

(file #21957)

Rik <rik5>
Group administrator
Sun 07 Nov 2010 06:51:21 PM UTC, original submission:  

Matlab and Octave 3.2.4 produce the following:
+verbatim
mod(2.1,0.1)
ans = 0
-verbatim-

The development version of Octave produces:

mod(2.1,0.1)
ans =  0.100000


This bug was initially reported under bug #31514.  It is being broken out because it is a different issue than that reported under #31514.

Rik <rik5>
Group administrator

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #21957:  fmodtst.cpp added by rik5 (601B - text/x-c++src)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by rik5 (Submitted the item)
  • -email is unavailable- added by rik5
  • -email is unavailable- added by rik5
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 5 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2010-11-09 jwe StatusNone Fixed
        Open/ClosedOpen Closed
    2010-11-07 rik5 Attached File- Added fmodtst.cpp, #21957
    2010-11-07 rik5 Carbon-Copy- Added kodo
        Carbon-Copy- Added jwe

    Back to the top

    Powered by Savane 3.13-caa5.
    Corresponding source code