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)

Comment Type & Canned Response:


Wed 08 Jul 2015 10:19:27 AM UTC, comment #4: 

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>
Project 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] )

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

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

(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 Caliari <caliari>
Project 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:

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 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 project members can vote.





    Follow 4 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    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.9