Hi @Hoibai,
This is a typical example of rounding errors due to numeric representation issues: Most numeric values with decimal places don't have an exact finite binary representation. In your example, the internal binary 64-bit floating-point representations of both 1604202898189.95 and 168579411743.91 are affected by rounding errors. Doing arithmetic with rounded values is likely to produce a rounding error in the result as well. As a consequence, a very small difference occurs between expressions that would be equal if the calculation was done in the decimal system (or with infinitely many bits in the binary system).
The same effect can occur with much smaller numbers:
11 data _null_;
12 if 1.86 ne 0.95 + 0.91 then put 'Surprised?';
13 run;
Surprised?
NOTE: DATA statement used (Total process time):
To avoid those surprises, use the ROUND function with an appropriate rounding unit.
Example:
14 data _null_;
15 if 1772782309933.86 = round(1604202898189.95 + 168579411743.91, 0.01) then put 'OK';
16 run;
OK
NOTE: DATA statement used (Total process time):
Here's how SAS computes the sum internally:
Looking at the internal binary representations (BINARY64. format) of 1604202898189.95 and 168579411743.91, we can write the sum mathematically in units of 2**40 in the binary system as follows:
1.0111010110000001111100011010001100001101111100110011001100110011001100110011001100...
+ 0.0010011101000000000111011010111100011111111010001111011111000010100011110101110000...
= 1.1001110011000010000011110101001000101101110111000010011(...)
The internal binary representations of numeric values (on Windows and Unix systems) use 52 mantissa bits plus one implied bit for the actual binary digits (while 11 bits are used for the exponent and 1 bit for the sign). In the first of the three lines above the 53 binary digits of 1604202898189.95 are shown in boldface: 41 of these, written in blue color, represent the integer part of the number, i.e., 1604202898189. The remaining 12 bits represent the fractional part (.95) -- more precisely: only a small portion of the infinitely many binary digits needed for an exact representation. The bits that don't fit into the total width of 53 bits are written in medium gray color. You can see the repeating pattern "1100" (underlined) of the periodic fraction 1/5=0.2 (binary: 0.00110011001100...), which together with the 0.75 (binary: 0.11) forms the 0.95. So this number has been rounded off internally.
The second line shows 168579411743.91 in a similar fashion. Note the three leading zeros due to the fact that the order of magnitude of this number is only 2**37 as opposed to 2**40. This allows for 53-38=15 bits representing the fractional part (.91), written in black (and red). But, again, even 15 bits are totally insufficient to accommodate the infinitely many bits required for an exact representation. I have added so many of the lost digits in gray that the (underlined) repeating pattern 10100011110101110000 (this time of length 20 rather than 4) can be recognized. Note, however, that the red "1" is rounded up from the original "0" (cf. the light green zero twenty digits later) because of the subsequent bits "11100...".
Now we can start the addition (like in elementary school), the result of which is shown in the third line: Remember that the gray digits are lost due to rounding, i.e., they are replaced by zeros altogether. So the first calculation adds the red "1" to the zero (which would be 1 without rounding) above it, resulting in the rightmost "1" in the third line, followed by another "0+1=1" for the bit left to it. The third-to-last bit is calculated as 0+0=0, the fourth-to-last bit as 1+1=10 (binary system!), resulting in the 0 (in the fourth-to-last bit of the third line) and a carryover of 1, which adds to the "1+1" coming next, and so on.
Finally, the result is rounded to fit into the 52+1 available bits. In the example above the last three bits "011" (gray) are rounded off. This matches the result of 1604202898189.95+168579411743.91, written in BINARY64. format, but not 1772782309933.86, written in BINARY64. format, which (in units of 2**40) would correspond to
1.1001110011000010000011110101001000101101110111000011
Obviously, the rounding errors coming from numeric representation in 64-bit floating-point format and from the calculation (rounding the result again) led to an incorrect least significant (i.e. rightmost) bit, representing 2**−12=0.000244140625 in this case. This is the "very small difference" mentioned earlier (well, "very small" compared to the large numbers involved), which causes the IF condition to be FALSE.
Applying the ROUND function with a rounding unit of 0.01 or 0.001, much greater than that small difference, corrects the error.
... View more
