[R] odd behaviour of %%?

[(Ted Harding)]

On 22-Nov-06 Jonathan Williams wrote:
> Dear R Helpers,
> 
> I am trying to extract the modulus from divisions by a
> sequence of fractions.
> I noticed that %% seems to behave inconsistently (to my
> untutored eye), thus:
> 
>> 0.1%%0.1
> [1] 0
>> 0.2%%0.1
> [1] 0
>> 0.3%%0.1
> [1] 0.1
>> 0.4%%0.1
> [1] 0
>> 0.5%%0.1
> [1] 0.1
>> 0.6%%0.1
> [1] 0.1
>> 0.7%%0.1
> [1] 0.1
>> 0.8%%0.1
> [1] 0
>> 0.9%%0.1

This is yet another manifestation of the fact that R (like
most computer languages) stores numbers as binary representations
of fixed length. This is OK for inetegers up to a certain
maximum value, but relatively few factions (those which are
multiples of powers of 1/2) will be stored exactly.

In particular, 0.1 is not stored exactly.

When you do 0.2 %% 0.1, you will be getting the remainder of
twice the binary representation of 0.1, modulo 0.1, which will
be zero; and similarly for 0.4 and 0.8, since these are going
up by powers of 2 and the relationships between their binary
representations will match what you would expect.

However, the binary representation of 0.3 does not correspond
exactly to 3 times the binary representation of 0.1:

  0.3 - 3*0.1
  [1] -5.551115e-17

so 3*0.1 is represented by a binary fraction slightly greater
than the binary representation of 0.3, so 0.1 only "goes into"
0.3 twice, with a remainder slightly less than  0.1 which rounds
to 0.1 when displayed as a result:

  > 0.3 %% 0.1
  [1] 0.1
  > 0.1 - (0.3 %% 0.1)
  [1] 2.775558e-17

Similarly for 0.5, 0.6, 0.7 and 0.9.

> The modulus for 0.1, 0.2, 0.4 and 0.8 is zero, as I'd expect.
> But, the modulus for 0.3, 0.6, 0.7 and 0.9 is 0.1 - which I
> did not expect. I can see no obvious rule that predicts whether
> x%%0.1 will give an answer of 0 or 0.1. I could find no
> explanation of the way that %% works in the R manuals. So,
> I have 3 questions:-
> 
> 1) Why is the modulus of 0.3%%0.1 (and 0.5%%0.1 and 0.6%%0.1...)
> not zero?

See above.

> 2) Are there any algorithms in R that use the %% operator with
> fractional divisors in this way,

I don't know -- though others will!

> and do they know about its apparently inconsisten behaviour?

People who write algorithms should be aware of such effects of
binary representation, so a well-written algorithm will take
account of this danger.

> 3) If %% is not intended for use with fractional divisors,
> then would it be a good idea to trap attempts to use them?

The best protection against this kind of thing is circumspection
on the part of the programmer -- who, therefore should do his
own trapping.

An alternative aproach to the calculations you made could be

  > 0.3 - 0.1*(0.3/0.1)
  [1] 0

which is "algebraically" equivalent, and gives the right answer:

  > identical(0.3 - 0.1*(0.3/0.1),0)
  [1] TRUE

However, this is not always going to work either:

  > 1000001.1 - (1000001.1/0.1)*0.1
  [1] -1.164153e-10

so it's prudent to wrap it in a suitable round():

  > round(1000001.1 - (1000001.1/0.1)*0.1, digits=9)
  [1] 0

but you still have to be aware of what value to set the number
of digits to round to:

  > round(1000001.1 - (1000001.1/0.1)*0.1, digits=10)
  [1] -1e-10

Hoping this helps,
Ted.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 22-Nov-06                                       Time: 11:33:05
------------------------------ XFMail ------------------------------