Rounding to multiples of a value lead to wrnog values, a small delta is added

[Georg]

I am using SAS 9.4.
The following code demonstates an error of the round (x, b) function:
Data work.tst;
uiva=-0.07918;
dciva=0.20066;
changeLines=( round(DCIVA,0.02) - round(UIVA,0.02))*10;
changeLines2=( round(DCIVA*50) - round(UIVA*50))/5;
changeLetters=changeLines*5;
if changelines=2.8 then a='Y'; else a='n';
if changelines2=2.8 then a2='Y'; else a2='n';
diff1=changelines - 2.8;
diff2=changelines2 - 2.8;
output;
run;

proc print data=Work.tst;
run;

Obs	uiva	dciva	changeLines	changeLines2	changeLetters	a	a2	diff1	diff2
1	-0.07918	0.20066	2.8	2.8	14	n	Y	4.4409E-16	0

So rounding by round(x, 0.02) leads to different result than round(x*50)/50.
The difference between both is a small delta=4.44e-16 !!!

This is a bad mistake, as typically it is often not seen when fomatted values are shown.
But it is effective when you make comparsions (see n for comparison a)

Is this bug going to be fixed somehow?

[Kurt_Bremser]

Differences like these are natural because of the way numeric data is stored and calculated. With 8-byte real, only about 15 decimal digits can be processed precisely. And most decimal fractions can not be converted to finite binary numbers anyway.

You will find the same "bug" in all software running on computers with a limited number of bits.

[ShenQicheng · Accepted solution]

I think the following document explains your question well.

SAS Help Center: ROUND Function

[Georg]

Thank you that makse it clear.

So always comaprind with delta:  e.g.   x1-x2 < 1e-10 
or maybe via round(x1, a) =round(x2,b)...

[Georg]

Let me add something:
The problem occurs only in the difference and not, if one calculates a single term:

Data work.tst;
uiva=-0.07918;
changeLines=round(UIVA,0.02);
changeLines2=round(UIVA*50)/50;
if changelines=-0.08 then a='Y'; else a='n';
if changelines2=-0.08 then a2='Y'; else a2='n';
diff1=changelines + 0.08;
diff2=changelines2 + 0.08;
output;
run;

Obs	uiva	changeLines	changeLines2	a	a2	diff1	diff2
1	-0.07918	-0.08	-0.08	Y	Y	0	0

That is quite strange! 

[FreelanceReinh]

Hello @Georg,

The ROUND function works well. As Kurt_Bremser said already, the issue is about arithmetic calculations with numbers of limited precision (due to numeric representation issues).

93    data _null_;
94    uiva=-0.07918;
95    dciva=0.20066;
96    if round(DCIVA,0.02)=0.2  then put 'OK 1';
97    if round(UIVA,0.02)=-0.08 then put 'OK 2';
98    if round(DCIVA*50)=10     then put 'OK 3';
99    if round(UIVA*50)=-4      then put 'OK 4';
100   if 0.2-(-0.08)=0.28       then put 'OK 5';
101   if 0.28*10 ne 2.8         then put 'Surprised?';
102   run;

OK 1
OK 2
OK 3
OK 4
OK 5
Surprised?

The first four OKs in the log above are what you would expect. The fifth OK is already a bit lucky because with slightly different numbers the comparison fails:

125   data _null_;
126   if 0.2-(-0.04) ne 0.24 then put 'Surprised again?';
127   run;

Surprised again?

This is due to the internal 64-bit floating-point binary representations of both 0.2 and 0.04 (and also 0.08) being inevitably affected by rounding errors as the exact binary representations of those numbers would require infinitely many binary digits -- which the 64 available bits obviously cannot accommodate. Arithmetic calculations with rounded numbers can easily propagate rounding errors and hence lead to slightly incorrect results. Sometimes, however, two rounding errors cancel out and the result is correct (as in 0.2-(-0.08)=0.28).

The "surprise" caused by the calculation 0.28*10 is just another example of a rounding error -- here contained in the internal representation of 0.28 only -- propagating to the result of a calculation. Let's take a closer look at this example and why it does not yield exactly 2.8.

The decimal number 0.28 written (mathematically) in the binary system is a repeating fraction:

0.01000111101011100001010001111010111000010100011110101110000101...

The 20-digit-pattern highlighted in blue is actually repeated forever (only three blocks of these are shown above). SAS uses the internal binary representation of 0.28 shown below, which can be seen by applying the BINARY64. format to 0.28.

0011111111010001111010111000010100011110101110000101000111101100

We recognize two complete copies of the repeating 20-digit-pattern, but the third copy is rounded (up) to fit into the 64 bits, more precisely: into the 52 mantissa bits plus the implied bit (highlighted in green above).

Multiplying this number by the decimal number 10 (binary: 1010) amounts to summing two shifted copies of it (because of the two 1s in "1010").

Written mathematically as multiples of 2**-2=0.25:

 1000.1111010111000010100011110101110000101000111101100  
+  10.001111010111000010100011110101110000101000111101100
---------------------------------------------------------
 1011.001100110011001100110011001100110011001100110011100

Now we see how the rounding error mentioned above affects the summation and leads to a result that is slightly too large. Moreover, it contains 55 binary digits, so it must be rounded again to 52+1=53 bits: The two trailing zeros at the end are cut off, but the rounding error persists in the remaining "1". The correct internal representation of 2.8 has a zero there, rounded down from (the repeating 4-digit pattern) 0011..., as can be seen in the BINARY64. format. The place value of the final incorrect "1" is 2**-51=4.44089...E-16, which is exactly the rounding difference that you found in variable diff1.

We can correct the rounding error by using the ROUND function with a suitable rounding unit such as 1E-9 (or anything else in {1E-2, 1E-3, ..., 1E-15}).

186   data _null_;
187   if round(0.28*10,1e-9)=2.8 then put 'OK';
188   run;

OK