[R] Best way to compute a sum

[Bernardo Rangel Tura]

On Sat, 2010-07-03 at 17:27 -0700, Roger Deangelis wrote:
> Hi Bernado,
> 
> In many financial applications if you convert the dollars and cents to
> pennies( ie $1.10 to 110) and divide by 100 at the vary end you can get
> maintain higher precision. This applies primarily to sums. This is similar
> to keeping track of the decimal fractions which have exact representations
> in floating point. It is good to know where the errors come from.

OK 
> 
> I suppose you could also improve the sum by understanding the decimal to
> binary rounding algorithms. I have noticed differences in computations
> between Sun hardware(FPUs) and Intel(FPUs). When presented the same set of
> operations, it appears that the floating point hardware do not do the
> operations in the same order. 

It's true
>  
> 
> Consider
> 
>    sprintf("%a", sum(rep(1.1,100)))     # overestimates the sum
> "0x1.b800000000001p+6"
>   "0x1.b800000000001p+6"
> 
>   sprintf("%a", sum(rep(11,100))/10)   # this gives the correct answer
> "0x1.b8p+6"
>   110 = 10 1110 = 32 + 8 + 4 + 2 - "0x1.b8p+6" - tricky due to 53 bit and
> little endian (I think this is right)
>   "0x1.b8p+6"
> 
>   note
>    cmp <- ifelse (sum(rep(1.1,100))==sum(rep(11,100))/10, "equal",
> "unequal")
>  
> [1] "unequal"
> 
>    cmp <- ifelse (sum(rep(1.1,100))>sum(rep(11,100))/10, "greater than",
> "less than orequal")
>   
> [1] "greater than"


Hi Roger,

First of all, is true 

> sum(rep(1.1,100))==sum(rep(11,100))/10
[1] FALSE

But is true too

> all.equal(sum(rep(1.1,100)),(sum(rep(11,100))/10))
[1] TRUE


I think you need read about "Guard Digits approach" 


-- 
Bernardo Rangel Tura, M.D,MPH,Ph.D
National Institute of Cardiology
Brazil