bugGNU Octave - Bugs: bug #45481, rem and fmod may give very wrong...

 
 

bug #45481: rem and fmod may give very wrong results for large arguments

Submitter:  Alexander Klein <matalex>
Submitted:  Mon 06 Jul 2015 10:23:50 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Need Info Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * dev Operating System:  * BSD
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Fri 12 Jan 2024 11:52:17 PM UTC, comment #5: 

Adding jwe to the CC list.

Octave produces nonsensical results for mod/rem operations on extremely large numbers.  Tracing this, I find that Octave has handrolled our own mod/rem functions in lo-mappers.h.  For example,


template <typename T>
T
mod (T x, T y)
{
  T retval;

  if (y == 0)
    retval = x;
  else
    {
      T q = x / y;

      if (x_nint (y) != y
          && (std::abs ((q - x_nint (q)) / x_nint (q))
              < std::numeric_limits<T>::epsilon ()))
        retval = 0;
      else
        {
          T n = std::floor (q);

          // Prevent use of extra precision.
          volatile T tmp = y * n;

          retval = x - tmp;
        }
    }

  if (x != y && y != 0)
    retval = copysign (retval, y);

  return retval;
}


Is there any reason to do this?  It seems better to rely on C/C++ standard libraries.  If I do, I find that


fmod (1e62, 10) = 4


which is not accurate in a pure math way (should be 0), but is consistent with Matlab and Sage and is not purely silly like 1.1e46.

Can we just map these functions onto the standard libraries as we do with trunc(), fix(), floor(), etc.?




Rik <rik5>
Group administrator
Wed 08 Jul 2015 10:19:27 AM UTC, comment #4: 

Alexander:
Right! When in doubt about floating point and many
other things lookup Kahan!

You may know that way back in the 1960's IBM started
shipping 360's with floating point without guard
digits. Kahan made them fix that at a very substantial
cost to IBM. It did not fix all the problems with 360 floating
point, but it helped a lot and was really the only
realistically possible improvement.

Michael Godfrey <godfrey>
Group Member
Wed 08 Jul 2015 09:22:19 AM UTC, comment #3: 

The problem is similar to #42627, but it goes in the opposite direction:

While the problem there is that rounding may have nudged a number away from being an integer, we have in our case numbers which are so large that their being or not being integers cannot be decided, anyway.

Look at the following:

octave:233> printf ( "%62.5f\n", 1e62 ./ [2 3 5 7 11 13 17 19 23] )
50000000000000001751099842971580586523040158899155912802435072.00000
33333333333333336403730152255666899093074731865430790378618880.00000
20000000000000001842238091353400139455844839119258474227171328.00000
14285714285714287437955120311274351342256590642183234843049984.00000
9090909090909090967973118230112855539046216263574686621761536.00000
7692307692307693235821987859916958415061241115021741838565376.00000
5882352941176471130070026868647099839954364446840727713873920.00000
5263157894736842064234558306593764730064337339464548234559488.00000
4347826086956522108589852626696554641307878733545426214977536.00000

Considering that 1e62 = 2^62 * 5^62, it is obvious that almost all of the above results are wrong in apparently being integers. Even long doubles in this region are so sparse, already, that integrity as a concept just crumbles away, but most of the time we can just ignore this fact, or are led astray by seeing decimal representations of binary encodings.

(Note: You may want to have a look at Floating-Point Arithmetic Besieged by “Business Decisions” and How Futile are Mindless Assessments of Roundoff in Floating-Point Computation ? by William Kahan, they are both quite amusing to read.)

Alexander Klein <matalex>
Wed 08 Jul 2015 08:23:29 AM UTC, comment #2: 

I'm reasonably confident that 4 is the correct answer for doubles.

The result by fmod_test has rounding errors creeping up above 1e22, but the remainders produced by fmod are consistent with the output of printf, and of lower magnitude than the divisor, at least.

I get 4 as the last digit with clang on FreeBSD as well as with gcc on OSX, but it may be worth noting that things are different for floats or long doubles, see attached program.

For OSX I get the following result, which is likely correct:

OSX% ./a.out
inf
inf
100000000000000003502199685943161173046080317798311825604870144
100000000000000003502199685943161173046080317798311825604870144
100000000000000000000982689773853339824972548992519069165944832
100000000000000000000982689773853339824972548992519069165944832

On FreeBSD, things go haywire, somewhat, insofar as the result of powf should have been Inf because it is larger than MAXFLOAT, and the result of powl is equal to the result of pow:

$ cc conv_test.c -lm
/tmp/conv_test-cca105.o: In function `main':
conv_test.c:(.text+0xf9): warning: powl has lower than advertised precision

$ ./a.out
inf
1000000030094932666179617348410047823344959136071346133401600
100000000000000003502199685943161173046080317798311825604870144
100000000000000003502199685943161173046080317798311825604870144
100000000000000000000982689773853339824972548992519069165944832
100000000000000003502199685943161173046080317798311825604870144

(file #34391)

Alexander Klein <matalex>
Wed 08 Jul 2015 07:02:35 AM UTC, comment #1: 

Hi, what we do expect here? The true result is 0. fmod_test gives 6. Matlab gives 4. Sage (mod(1e62,10)) gives 4.

Isn't it similar to #42627? If x / y is assumed integer (for any y, not only if y is not integer), the result should be 0.

Marco

Marco Caliari <caliari>
Group Member
Mon 06 Jul 2015 10:23:50 AM UTC, original submission:  

Octave 3.6.4 on OSX and 3.8.2 on FreeBSD may give completely wrong results when using rem/fmod:


octave:191> rem(1e62,10)
ans =  1.1418e+46


This seems to be caused by the miracles of floating point in conjunction with literally using the formula from the docs:


octave:196> x=1e62, y=10, x - y .* fix (x ./ y)
x =  1.0000e+62
y =  10
ans = -1.1418e+46


The correct result might be computed much faster using the fmod-function from libm - see attached program - but maybe there was a reason not to use it, in the first place?

Alexander Klein <matalex>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #34391:  conv_test.c added by matalex (331B - application/octet-stream)
file #34384:  fmod_test.c added by matalex (248B - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by rik5
  • -email is unavailable- added by mtmiller (Updated the item)
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by matalex (Submitted the item)
  •  

    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
    2024-01-12 rik5 Carbon-Copy- Added jwe
    2018-05-01 mtmiller StatusNone Need Info
        Release3.8.2 dev
    2015-07-08 matalex Attached File- Added conv_test.c, #34391
    2015-07-06 matalex Attached File- Added fmod_test.c, #34384

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code