Hello!
I understand a little bit about floating point storage; however, could someone explain these results? I can understand a number not being stored exactly due the the precision of the system and how that leads to conflict in calculations sometimes. I'm a little confused how these two representations of the same number that do not include any calculations get stored differently. Is SAS turning scientific notation into a calculation rather than just moving the decimal place over?
Warm regards,
Michael
3452 data _null_;
3453 if 727 = 7.27E2 then put 'Equal'; else put 'Unequal';
3454 if 72.7 = 7.27E1 then put 'Equal'; else put 'Unequal';
3455 if 7.27 = 7.27E0 then put 'Equal'; else put 'Unequal';
3456 if 0.727 = 7.27E-1 then put 'Equal'; else put 'Unequal';
3457 if 0.0727 = 7.27E-2 then put 'Equal'; else put 'Unequal';
3458 if 0.00727 = 7.27E-3 then put 'Equal'; else put 'Unequal';
3459 if 0.000727 = 7.27E-4 then put 'Equal'; else put 'Unequal';
3460 if 0.0000727 = 7.27E-5 then put 'Equal'; else put 'Unequal';
3461 put '---';
3462 if 0.0000727 = round(7.27E-5 , 1E-16) then put 'Equal'; else put 'Unequal';
3463 if 0.0000727 = round(0.0000727, 1E-16) then put 'Equal'; else put 'Unequal';
3464 if 7.27E-5 = round(7.27E-5 , 1E-16) then put 'Equal'; else put 'Unequal';
3465 run;
Equal
Equal
Unequal
Equal
Equal
Equal
Equal
Unequal
---
Equal
Equal
Unequal
Here's another example that I don't understand.
3701 data _null_;
3702 if 7900000000000000000 = round(7900000000000000000, 1000000000 ) then put 'Equal'; else put 'Unequal';
3703 if 7900000000000000000 = round(7900000000000000000, 10000000000 ) then put 'Equal'; else put 'Unequal';
3711 run;
Equal
Unequal
@FreelanceReinh wrote:
...
The results of the DO loops in @mkeintz's code can be explained as follows:
Positive integer multiples of 1E9 (1E10) with up to 62 (63) digits in their mathematical binary representation have an exact internal representation in SAS under Windows. This is because the >=9 (10) trailing binary zeros can safely be truncated when the 61-bit (62-bit) mantissa is shortened to the 52 available bits. In the binary system 7.9E18 has 63 digits. So, the second DO loop operates entirely in the "safe zone" where no rounding errors occur.
The course of events in the first DO loop is different: After 4,611,686,019 iterations, i.e., more than 3 billion iterations before approx. 7.9E18 is reached, xnew exceeds 2**62=4.61...E18 and the first of millions (if not billions) of rounding errors occurs: 4611686019000000512. These rounding errors accumulate dramatically, leading to a substantial deviation from 7.9E18 in the final value of xnew. Moreover, the loop does not terminate after 7,900,000,000, but only after 7,900,001,684 iterations.
More generally on the subject of precision for numbers > CONSTANT('exactint',8) (call it EI8, which equals =9,007,199,254,740,992).
While double-precision floating point can exactly represent every positive integer <=EI8, it can only exactly represent this subset of integers above EI8:
BTW, this goes downward too, so that you have exact representation of these values below EI8:
Of course EI8 is a different value for SAS on ibm mainframe's, which I recall use fewer bits for the exponent and more for the mantissa.
Here's the promised program log. It forces the SAS compiler to assign floating point values for literals from 9,007,199,254,740,992 through 9,007,199,254,741,008. So all odd numbers are set to even numbers, but notice half the odd numbers are rounded up, half are rounded down, in a regular pattern. I believe this is some IEEE standard meant to avoid bias in "rounding". The result is that half of the even number get no rounding results, and the other half get all of the rounding results.
2390 data _null_;
2391 x0992=9007199254740992;
2392 x0993=9007199254740993;
2393 x0994=9007199254740994;
2394 x0995=9007199254740995;
2395 x0996=9007199254740996;
2396 x0997=9007199254740997;
2397 x0998=9007199254740998;
2398 x0999=9007199254740999;
2399 x1000=9007199254741000;
2400 x1001=9007199254741001;
2401 x1002=9007199254741002;
2402 x1003=9007199254741003;
2403 x1004=9007199254741004;
2404 x1005=9007199254741005;
2405 x1006=9007199254741006;
2406 x1007=9007199254741007;
2407 x1008=9007199254741008;
2408 put (_all_) (=comma21.0 /);
2409 run;
x0992=9,007,199,254,740,992
x0993=9,007,199,254,740,992
x0994=9,007,199,254,740,994
x0995=9,007,199,254,740,996
x0996=9,007,199,254,740,996
x0997=9,007,199,254,740,996
x0998=9,007,199,254,740,998
x0999=9,007,199,254,741,000
x1000=9,007,199,254,741,000
x1001=9,007,199,254,741,000
x1002=9,007,199,254,741,002
x1003=9,007,199,254,741,004
x1004=9,007,199,254,741,004
x1005=9,007,199,254,741,004
x1006=9,007,199,254,741,006
x1007=9,007,199,254,741,008
x1008=9,007,199,254,741,008
As to your first set of examples: I hadn't thought about it before, but I guess I'm not surprised that SAS (and other compilers?) might interpret real value literals differently when specified in scientific notation - even when the intended values are identical. I think practitioners with other compilers have best practices when it comes to specifying numeric literals stored as floating point.
But I suspect you'll have to ask SAS Institute about that one.
As to the second example, there must be something about how the compiler sends the division (multiplication by inverse I suppose) request to the cpu, because the round functions produce the R1 and R2 values below:
1124 data _null_;
1125 e8=constant('exactint',8);
1126 r1=round(7900000000000000000, 1000000000 );
1127 r2=round(7900000000000000000, 10000000000 );
1128 put (e8 r1 r2) (=21.0 / );
1129 run;
e8=9007199254740992
r1=7900000000000000000
r2=7899999999999998976
R2 is exactly 1,024 less than R1. This makes sense because none of the integers between R1 and R2 are representable in floating point double-precision (8 byte) storage (note the value for E8, which is the largest consecutive integer representable in floating point double precision).
I guess the division which may be the first stage of the round function generates different results.
I did this program, which starts from zero and adds 1000000000 (and in a second program 10000000000) at a time, to print the nearest values to =7900000000000000000. Notice the first sum never exactly matches 7900000000000000000:
1307 data _null_;
1308 xnew=0;
1309 do until (xnew>=7900000000000000000);
1310 xold=xnew;
1311 xnew+1000000000;
1312 end;
1313 put (xold xnew) (=21.0 / );
1314 run;
xold=7899999999382380544
xnew=7900000000382380032
NOTE: DATA statement used (Total process time):
real time 16.72 seconds
cpu time 16.92 seconds
1315
1316 data _null_;
1317 xnew=0;
1318 do until (xnew>=7900000000000000000);
1319 xold=xnew;
1320 xnew+10000000000;
1321 end;
1322 put (xold xnew) (=21.0 / );
1323 run;
xold=7899999990000000000
xnew=7900000000000000000
NOTE: DATA statement used (Total process time):
real time 1.66 seconds
cpu time 1.65 seconds
Hello @Kastchei,
Many thanks for providing these examples. I will add the 7.27 ne 7.27E0 case to my collection of numeric representation oddities. One of my standard examples in this regard is 0.00001 ne 1.0E-5. In 2016 there was a long discussion (involving relevant SAS personnel) about cases where adding insignificant trailing zeros to a decimal fraction (like 0.500000...) changes the internal floating point representation: https://communities.sas.com/t5/SAS-Programming/Expected-numeric-precision-behaviour-or-unexpected-is... (which led to a supplement of the SAS documentation). I think for a satisfactory explanation of both types of examples one would need to know and understand the (assembly) routines that SAS uses internally for the "decimal-to-binary" conversion. (The calculation 7.27*1E0 cannot explain the excessive rounding error contained in 7.27E0.)
I'll take a closer look at the examples involving the ROUND function tomorrow (or rather later today, it's after midnight here) and then continue this post. It's possible that what is known about the internal workings of the ROUND function can explain some of these.
[Hours later ...]
The internal binary floating-point representation of 0.0000727 is what I would expect by IEEE-754 standards. The least significant bit has a place value of 2**-66 = 1.355E-20. The numeric representation error (incurred by rounding up) amounts to 5.286E-21. The internal representation of 7.27E-5, however, is "incorrectly" rounded down, leading to a larger representation error (but with negative sign): -8.266E-21. Rounding either of the two numbers to the closest multiple of 1E-16 (a rounding unit much larger than 1.355E-20) results in 0.0000727, as it should. This explains (or: is consistent with) what you obtained from the code in lines 3462-3464 of your log.
The case of 7900000000000000000 (=7.9E18) shown in lines 3702-3703 of your log is more mysterious to me. First of all, the internal representation of this number is exact, despite its size (beyond constant('exactint')), i.e., representation error is zero. This is because the 17 decimal zeros (factor 10**17 = 5**17 * 2**17) entail 17 binary zeros (factor 2**17) in the mathematical binary representation, which reduces the number of significant mantissa bits from 62 to 45 (excl. the implied bit), while 52 bits are available (under Windows). Also, SAS agrees that 7900000000000000000 = 7900000000*1E9 = 790000000*1E10. Still, the result of the ROUND function using rounding unit 1E10 is incorrect (too small) in the least significant bit, which represents 2**(62-52)=2**10=1024.
Unfortunately, the ROUND function documentation provides only hints about the internal workings of the function. It suggests that 0.5 times the rounding unit plus a small "epsilon" is added to the first argument and the result is then truncated (roughly speaking), which makes sense. That addition might already contribute to the error in the final result: Note that, unlike 7.9E18 itself, 7.9E18 + 0.5E10 = 7900000005000000000 (disregarding the "epsilon") does not have an exact internal representation: The representation error is -512. Similarly, the representation of 7.9E18 + 0.5E9 is not exact either, but the error happens to be smaller (in absolute value): -256.
However, I wasn't able to reproduce the incorrect result by means of data step calculations as suggested in the documentation:
366 data test; 367 do eps=0 to 1e7; 368 x=7900000005000000000+eps; 369 m=modz(x,1e10); 370 r=x-m; 371 if r ne 7900000000000000000 then output; 372 end; 373 run; NOTE: The data set WORK.TEST has 0 observations and 4 variables.
[EDIT: Addendum: I've just found out that the ROUNDZ function does not show that surprising behavior, which suggests that perhaps the fuzzing part of the ROUND algorithm causes the deviation from the expected result in this case.]
The results of the DO loops in @mkeintz's code can be explained as follows:
Positive integer multiples of 1E9 (1E10) with up to 62 (63) digits in their mathematical binary representation have an exact internal representation in SAS under Windows. This is because the >=9 (10) trailing binary zeros can safely be truncated when the 61-bit (62-bit) mantissa is shortened to the 52 available bits. In the binary system 7.9E18 has 63 digits. So, the second DO loop operates entirely in the "safe zone" where no rounding errors occur.
The course of events in the first DO loop is different: After 4,611,686,019 iterations, i.e., more than 3 billion iterations before approx. 7.9E18 is reached, xnew exceeds 2**62=4.61...E18 and the first of millions (if not billions) of rounding errors occurs: 4611686019000000512. These rounding errors accumulate dramatically, leading to a substantial deviation from 7.9E18 in the final value of xnew. Moreover, the loop does not terminate after 7,900,000,000, but only after 7,900,001,684 iterations.
@FreelanceReinh wrote:
...
The results of the DO loops in @mkeintz's code can be explained as follows:
Positive integer multiples of 1E9 (1E10) with up to 62 (63) digits in their mathematical binary representation have an exact internal representation in SAS under Windows. This is because the >=9 (10) trailing binary zeros can safely be truncated when the 61-bit (62-bit) mantissa is shortened to the 52 available bits. In the binary system 7.9E18 has 63 digits. So, the second DO loop operates entirely in the "safe zone" where no rounding errors occur.
The course of events in the first DO loop is different: After 4,611,686,019 iterations, i.e., more than 3 billion iterations before approx. 7.9E18 is reached, xnew exceeds 2**62=4.61...E18 and the first of millions (if not billions) of rounding errors occurs: 4611686019000000512. These rounding errors accumulate dramatically, leading to a substantial deviation from 7.9E18 in the final value of xnew. Moreover, the loop does not terminate after 7,900,000,000, but only after 7,900,001,684 iterations.
More generally on the subject of precision for numbers > CONSTANT('exactint',8) (call it EI8, which equals =9,007,199,254,740,992).
While double-precision floating point can exactly represent every positive integer <=EI8, it can only exactly represent this subset of integers above EI8:
BTW, this goes downward too, so that you have exact representation of these values below EI8:
Of course EI8 is a different value for SAS on ibm mainframe's, which I recall use fewer bits for the exponent and more for the mantissa.
Here's the promised program log. It forces the SAS compiler to assign floating point values for literals from 9,007,199,254,740,992 through 9,007,199,254,741,008. So all odd numbers are set to even numbers, but notice half the odd numbers are rounded up, half are rounded down, in a regular pattern. I believe this is some IEEE standard meant to avoid bias in "rounding". The result is that half of the even number get no rounding results, and the other half get all of the rounding results.
2390 data _null_;
2391 x0992=9007199254740992;
2392 x0993=9007199254740993;
2393 x0994=9007199254740994;
2394 x0995=9007199254740995;
2395 x0996=9007199254740996;
2396 x0997=9007199254740997;
2397 x0998=9007199254740998;
2398 x0999=9007199254740999;
2399 x1000=9007199254741000;
2400 x1001=9007199254741001;
2401 x1002=9007199254741002;
2402 x1003=9007199254741003;
2403 x1004=9007199254741004;
2404 x1005=9007199254741005;
2405 x1006=9007199254741006;
2406 x1007=9007199254741007;
2407 x1008=9007199254741008;
2408 put (_all_) (=comma21.0 /);
2409 run;
x0992=9,007,199,254,740,992
x0993=9,007,199,254,740,992
x0994=9,007,199,254,740,994
x0995=9,007,199,254,740,996
x0996=9,007,199,254,740,996
x0997=9,007,199,254,740,996
x0998=9,007,199,254,740,998
x0999=9,007,199,254,741,000
x1000=9,007,199,254,741,000
x1001=9,007,199,254,741,000
x1002=9,007,199,254,741,002
x1003=9,007,199,254,741,004
x1004=9,007,199,254,741,004
x1005=9,007,199,254,741,004
x1006=9,007,199,254,741,006
x1007=9,007,199,254,741,008
x1008=9,007,199,254,741,008
Are you ready for the spotlight? We're accepting content ideas for SAS Innovate 2025 to be held May 6-9 in Orlando, FL. The call is open until September 16. Read more here about why you should contribute and what is in it for you!
Learn how use the CAT functions in SAS to join values from multiple variables into a single value.
Find more tutorials on the SAS Users YouTube channel.
Select SAS Training centers are offering in-person courses. View upcoming courses for: