Re: Accuracy, format and round

AIO

I understand that SAS performs arithmetic in a binary system, but I still have following questions about it:

Variable X1 X2: since variable flag="Y", x1 should be less than x2. However, when variable x1 and x2 are displayed as best32. format, it's "73.75". How could I find the difference directly except using HEXw. format? Is there a way that let X1 display as "73.7499999999999999999999999999999999"?
Variable Z1 Z2: actually, x1 should be less than 73.75. It's specified that "When the value to be rounded is approximately halfway between two multiples of the rounding unit, the ROUND function rounds up the absolute value and restores the original sign" in ROUND function. In this case, it should be 73.7 when I use ROUND function to keep 1 decimal place. However, Z1 is 73.8, why?
⁠Variable W1 W2: I want to know what's rule when using PUT function to keep the decimal place? Is it "When the value to be rounded is approximately halfway between two multiples of the rounding unit, the ROUND function rounds up the absolute value and restores the original sign"?

data test;
	x1=(36.9-31)*12.5;
	x2=73.75;
	if x1 <x2 then flag="Y";

	y1=x1;
	y2=x2;
	format x1 x2 best32. y1 y2 hex16.;

	z1=round(x1,0.1);
	z2=round(x2,0.1);

	w1=strip(put(x1,5.1));
	w2=strip(put(x2,5.1));
run;

Ksharp

1)

For your first question, sas is unable to store 73.75 exactly in computer sometimes, so you can't always expect the expression is equal. If you want to compare them ,the best way is using ROUND() function.

if you need to display "73.7499999999999999999999999999999999" , you need to resort to PICTURE format.

proc format;
picture fmt
low-high='09.99999999999999'
;
run;

data test;
	x1=(36.9-31)*12.5;
	x2=73.75;
	if x1 <x2 then flag="Y";
put x1=  fmt32.20 x2= fmt32.20  flag=;

    call missing(flag);
	if round(x1,0.1) <round(x2,0.1) then flag="Y";
put x1=  32.20 x2= 32.20 flag=;



    z1=round(x1,0.1);
	z2=round(x2,0.1);
put z1=  32.20 z2= 32.20;

    w1=strip(put(x1,5.1));
	w2=strip(put(x2,5.1));
put w1= w2=;
run;

x1=73.74999999999999 x2=73.75000000000000 flag=Y
x1=73.75000000000000000000 x2=73.75000000000000000000 flag=
z1=73.80000000000000000000 z2=73.80000000000000000000
w1=73.7 w2=73.8
NOTE: 数据集 WORK.TEST 有 1 个观测和 7 个变量。
NOTE: “DATA 语句”所用时间（总处理时间）:
      实际时间          0.03 秒
      CPU 时间          0.03 秒

2)

For your second question, in sas there is a constant epsolon:

406  data _null_;
407  ep=constant('MACEPS');
408  put ep=;
409  run;

ep=2.220446E-16

Once the difference of two numbers is less than epsolon, then sas would take these two number are same, your situation is this.

3)

For your third question, I also get surprise at this result. I also think

put(x1,5.1)

would be like this

put(round(x1,5.1),5.1)

and would expect both are 73.8 . Not one is 73.7 and another is 73.8 .

I suggest you to contact with SAS -Technical -Support to discover what is wrong with PUT() .

AIO

For Q2, the difference between X1 and X2 is not less than machine epsilon. SAS shouldn't take them as same.

data test;
	x1=(36.9-31)*12.5;
	x2=73.75;
	if x1 <x2 then flag="Y";
	diff=x2-x1;
	ep=constant('MACEPS');
	if diff<ep then lsfl="Y";
run;

FreelanceReinh

@AIO wrote:

For Q2, the difference between X1 and X2 is not less than machine epsilon. SAS shouldn't take them as same.

See the paragraph in the documentation of the ROUND function starting with "The approximation is relative to the size of the value to be rounded ..." The "fuzzing" of the ROUND function does not use an absolute additive constant (such as machine epsilon).

FreelanceReinh

Hello @AIO,

I agree with @Ksharp that applying the ROUND function to the result of calculations like that in variable X1 would be best practice as it avoids propagating the unavoidable rounding error (due to numeric representation issues).

While 73.75 does have an exact, finite representation in the binary system (1001001.11), your rounding issue comes from the calculation: The result of 36.9-31, computed in the binary system, using the internal binary floating-point representations of 36.9 (which already contains a numeric representation error) and 31, is slightly less than 5.9.

170   data _null_;
171   d=5.9-(36.9-31);
172   put d;
173   run;

1.776357E-15

This difference, 2**-49=1.776356839...E-15, reflects a binary digit, actually three digits 010, which are missing (compared to the internal representation of 5.9) at the end of the mantissa because the internal representation of 36.9 contains five rather than two bits for the integer part at the beginning of the mantissa. Of course, multiplying the result by 12.5 doesn't bring those lost bits back, hence the result <73.75.

To obtain more decimal places than the picture format suggested by @Ksharp, you could apply a user-defined function such as BASECONVFP and translate the internal binary representation of the fractional part of X1 to the decimal system:

183   data _null_;
184   x1=(36.9-31)*12.5;
185   b=put(x1, binary64.);
186   length dec $50;
187   dec=baseconvfp(substr(b,19),2,10,50);
188   put dec;
189   run;

74999999999998578914528479799628257751464843750000

So, adding the integer part, you get 73.7499999999999857891452847979962825775146484375.

As to your second question: To avoid the "fuzzing" of the ROUND function, use the ROUNDZ function:

z1=roundz(x1,0.1);
z2=roundz(x2,0.1);

Result: z1=73.7, z2=73.8.

Regarding your third question, I think the ROUND function uses a more sophisticated algorithm than the PUT function and hence is the preferred tool for rounding numeric values.

Accuracy, format and round