Scientific notation and numeric precision

[Kastchei]

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

[mkeintz]

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

[FreelanceReinh]

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.

[mkeintz · Accepted solution]

---

@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:

• Every even number up to 2*EI8.   All other values in this range will be rounded up or down to an even number (see the program log below).
• Every 0mod4 number up to 4*EI8
• Every 0mod8 number up to 8*EI8
• ...
• Every 0mod1024 number up to 1024*EI8 = EI8*(2**10), which happens to include 7,900,000,000.
• ...
• And I suspect every 0mod(2**1024) number up to CONSTANT('BIG') =1.797693134862300E+308

BTW, this goes downward too, so that you have exact representation of these values below EI8:

• Every integer in [EI8/2, EI8)
• Every integer plus every exact-half in [EI8/4, EI8/2)
• Every integer, exact half, plus every exact quarter in [EI8/8, EI8/4)   
• etc., etc.

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

[Kastchei]

Thank you both for the detailed explanations and the simulations! This was way more helpful than I was expecting.

I am comparing two ways to round to significant figures. One was extremely simple:

input(put(x, E9.), E9.)

The other used some log10 and other math to determine which decimal place is needed for the second argument of the round function.

I wasn't so surprised to see differences for extremely small values, but some of the ones in my original post, like 7.27E0, where surprising.

From my own experimenting, it seems that the scientific notation version works a little better for very large numbers: that is, it's exact for slightly more magnitudes of order, whereas the round seems to have a mixup along the way.

Thanks again. Happy holidays!

Warm regards,
Michael