[R] Rounding problem R vs Excel

[Marc Schwartz]

Hi all,

I thought that I would follow up on this thread with some "refined"
information.

As I mentioned in a prior post in this thread, I posted a query to the
OOo folks regarding the inability to replicate the IEEE representation
issue in OOo Calc.

Recall the following results:

OOo Calc 1.0.2 and 1.1 Beta2:

Cell Formula          Value
= 4.145 * 100 + 0.5   4.15000000000000000000E+02
= 0.5 - 0.4 - 0.1     0.00000000000000000000E+00
=(0.5 - 0.4 - 0.1)    0.00000000000000000000E+00


MS Excel 2002 (XP)

Cell Formula          Value
= 4.145 * 100 + 0.5   4.15000000000000000000E+02
= 0.5 - 0.4 - 0.1     0.00000000000000000000E+00
=(0.5 - 0.4 - 0.1)    -2.77555756156289000000E-17


Gnumeric 1.0.12:

Cell Formula          Value
= 4.145 * 100 + 0.5   +4.14999999999999943157E+02
= 0.5 - 0.4 - 0.1     -2.77555756156289135106E-17
*Gnumeric does not appear to allow the surrounding parens.


R 1.7.1 Beta:

> print(4.145 * 100 + 0.5, digits = 20)
[1] 414.99999999999994
> formatC(4.145 * 100 + 0.5, format = "E", digits = 20)
[1] "4.14999999999999943157E+02"

> print(0.5 - 0.4 - 0.1, digits = 20)
[1] -2.775557561562891e-17
> formatC(0.5 - 0.4 - 0.1, format = "E", digits = 20)
[1] "-2.77555756156289135106E-17"


Today, I received a reply from Niklas Nebel at Sun to indicate that
indeed OOo has included the same (or perhaps more correctly, a similar)
"optimization" that MS incorporated into Excel starting with Excel 97,
which renders results "close to zero" as zero.

In OOo Calc 1.1 Beta, Niklas indicated that the relevant source code is
in:

"sal/inc/rtl/math.hxx"

Since the full source tarball is rather large, a subset of the relevant
source code from that file is below. Note the use of the approxEqual()
function as the basis for subsequent arithmetic. 

This code should, in theory, offer some insight into what Excel is doing
as well, with the exception of course of the formula parsing issue,
which is clearly a bug.

HTH,

Marc Schwartz


***************

Copyright: 2002 by Sun Microsystems, Inc.

/** Test equality of two values with an accuracy of the magnitude of the
    given values scaled by 2^-48 (4 bits roundoff stripped).

    @ATTENTION
    approxEqual( value!=0.0, 0.0 ) _never_ yields true.
*/
inline bool approxEqual(double a, double b)
{
  if ( a == b )
    return true;
  double x = a - b;
    return (x < 0.0 ? -x : x)
            < ((a < 0.0 ? -a : a) * (1.0 / (16777216.0 * 16777216.0)));
}

/** Add two values.

If signs differ and the absolute values are equal according to
approxEqual() the method returns 0.0 instead of calculating the sum.

If you wanted to sum up multiple values it would be convenient not to
call approxAdd() for each value but instead remember the first value not
equal to 0.0, add all other values using normal + operator, and with the
result and the remembered value call approxAdd().
*/
inline double approxAdd(double a, double b)
{
  if ( ((a < 0.0 && b > 0.0) || (b < 0.0 && a > 0.0))
     && approxEqual( a, -b ) )
  return 0.0;
  return a + b;
}

/** Substract two values (a-b).

If signs are identical and the values are equal according to
approxEqual() the method returns 0.0 instead of calculating the
substraction.
*/
inline double approxSub(double a, double b)
{
if ( ((a < 0.0 && b < 0.0) || (a > 0.0 && b > 0.0)) && 
     approxEqual(a, b) )
  return 0.0;
  return a - b;
}

/** floor() method taking approxEqual() into account.

  Use for expected integer values being calculated by double functions.

    @ATTENTION
    The threshhold value is 3.55271e-015
*/

inline double approxFloor(double a)
{
  double b = floor( a );
  // The second approxEqual() is necessary for values that are near the
limit
  // of numbers representable with 4 bits stripped off. (#i12446#)
  if ( approxEqual( a - 1.0, b ) && !approxEqual( a, b ) )
    return b + 1.0;
  return b;
}

/** ceil() method taking approxEqual() into account.

Use for expected integer values being calculated by double functions.

    @ATTENTION
    The threshhold value is 3.55271e-015
*/
inline double approxCeil(double a)
{
  double b = ceil( a );
  // The second approxEqual() is necessary for values that are near the 
limit
  // of numbers representable with 4 bits stripped off. (#i12446#)
  if ( approxEqual( a + 1.0, b ) && !approxEqual( a, b ) )
    return b - 1.0;
  return b;
}