Add a New Comment (Rich Markup)

Comment Type & Canned Response: None None > Multiple Canned Responses Fixed in development Crash with no stack trace Already fixed in newer version Fixed in stable Bad description Bad description and crash Bad stack trace Obsolete version Duplicate and not fixed Duplicate and needs more info Duplicate and fixed Need info and old

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.

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.)

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)

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

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?

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

Attach Files: Comment:

Depends on the following items: None found

Items that depend on this one: None found

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.

Please enter the title of George Orwell's famous dystopian book (it's a date):

Follow 4 latest changes.

Copyright © 2022 Free Software Foundation, Inc. Verbatim copying and distribution of this entire article is permitted in any medium, provided this notice is preserved. The Levitating, Meditating, Flute-playing Gnu logo is a GNU GPL'ed image provided by the Nevrax Design Team. Source Code

Powered by Savane 3.9