Season's Greetings (and great news ... )!

I wouldn't blame R for floating-point arithmetic and our personal
feeling of what 'zero' should be.

options(digits=20)
pi

[1] 3.141592653589793116

sqrt(pi)^2

[1] 3.1415926535897926719

(pi - sqrt(pi)^2) < 1e-15

[1] TRUE

There was a similar post before, for example see:
http://r.789695.n4.nabble.com/Why-does-sin-pi-not-return-0-td4676963.html

There is an example by Martin Maechler (author of Rmpfr) on how to use
arbitrary precision
with your arithmetic.

On 22 December 2013 10:59, Ted Harding <Ted.Harding at wlandres.net> wrote:

Greetings All!
With the Festive Season fast approaching, I bring you joy
with the news (which you will surely wish to celebrate)
that R cannot do arithmetic!

Usually, this is manifest in a trivial way when users report
puzzlement that, for instance,

  sqrt(pi)^2 == pi
  # [1] FALSE

which is the result of a (generally trivial) rounding or
truncation error:

  sqrt(pi)^2 - pi
  [1] -4.440892e-16

But for some very simple calculations R goes off its head.

I had originally posted this example some years ago, but I
have since generalised it, and the generalisation is even
more entertaining than the original.

The Original:
Consider a sequence generated by the recurrence relation

  x[n+1] = 2*x[n] if 0 <= x[n] <= 1/2
  x[n+1] = 2*(1 - x[n]) if 1/2 < x[n] <= 1

(for 0 <= x[n] <= 1).

This has equilibrium points (x[n+1] = x[n]) at x[n] = 0
and at x[n] = 2/3:

  2/3 -> 2*(1 - 2/3) = 2/3

It also has periodic points, e.g.

  2/5 -> 4/5 -> 2/5 (period 2)
  2/9 -> 4/9 -> 8/9 -> 2/9 (period 3)

The recurrence relation can be implemented as the R function

  nextx <- function(x){
    if( (0<=x)&(x<=1/2) ) {x <- 2*x} else {x <- 2*(1 - x)}
  }

Now have a look at what happens when we start at the equilibrium
point x = 2/3:

  N <- 1 ; x <- 2/3
  while(x > 0){
    cat(sprintf("%i: %.9f\n",N,x))
    x <- nextx(x) ; N <- N+1
  }
  cat(sprintf("%i: %.9f\n",N,x))

Run that, and you will see that successive values of x collapse
towards zero. Things look fine to start with:

  1: 0.666666667
  2: 0.666666667
  3: 0.666666667
  4: 0.666666667
  5: 0.666666667
  ...

but, later on,

  24: 0.666666667
  25: 0.666666666
  26: 0.666666668
  27: 0.666666664
  28: 0.666666672
  ...

  46: 0.667968750
  47: 0.664062500
  48: 0.671875000
  49: 0.656250000
  50: 0.687500000
  51: 0.625000000
  52: 0.750000000
  53: 0.500000000
  54: 1.000000000
  55: 0.000000000

What is happening is that, each time R multiplies by 2, the binary
representation is shifted up by one and a zero bit is introduced
at the bottom end. To illustrate this, do the calculation in
7-bit arithmetic where 2/3 = 0.1010101, so:

0.1010101  x[1], >1/2 so subtract from 1 = 1.0000000 -> 0.0101011,
and then multiply by 2 to get x[2] = 0.1010110. Hence

0.1010101  x[1] -> 2*(1 - 0.1010101) = 2*0.0101011 ->
0.1010110  x[2] -> 2*(1 - 0.1010110) = 2*0.0101010 ->
0.1010100  x[3] -> 2*(1 - 0.1010100) = 2*0.0101100 ->
0.1011000  x[4] -> 2*(1 - 0.1011000) = 2*0.0101000 ->
0.1010000  x[5] -> 2*(1 - 0.1010000) = 2*0.0110000 ->
0.1100000  x[6] -> 2*(1 - 0.1100000) = 2*0.0100000 ->
0.1000000  x[7] -> 2*0.1000000 = 1.0000000 ->
1.0000000  x[8] -> 2*(1 - 1.0000000) = 2*0 ->
0.0000000  x[9] and the end of the line.

The final index of x[i] is i=9, 2 more than the number of binary
places (7) in this arithmetic, since 8 successive zeros have to
be introduced. It is the same with the real R calculation since
this is working to .Machine$double.digits = 53 binary places;
it just takes longer (we reach 0 at x[55])! The above collapse
to 0 occurs for any starting value in this simple example (except
for multiples of 1/(2^k), when it works properly).

Generalisation:
This is basically the same, except that everything is multiplied
by a scale factor S, so instead of being on the interval [0,1].
it is on [0,S], and

  x[n+1] = 2*x[n] if 0 <= x[n] <= S/2
  x[n+1] = 2*(S - x[n]) if S/2 < x[n] <= S
(for 0 <= x[n] <= S).

Again, x[n] = 2*S/3 is an equilibrium point. 2*S/3 > S/2, so

  x[n] -> 2*(S - 2*S/3) = 2*(S/3) = 2*S/3

Functions to implement this:

  nxtS <- function(x,S){
    if((x >= 0)&(x <= S/2)){ x<- 2*x } else {x <- 2*(S-x)}
  }

  S <- 6 ##  Or some other value of S
  Nits <- 100
  x <- 2*S/3
  N <- 1 ; print(c(N,x))
  while(x>0){
  if(N > Nits) break   ### to stop infinite looping
  N <- (N+1) ; x <- nxtS(x,S)
  print(c(N,x))
}

The behaviour of the sequence now depends on the value of S.

If S is a multiple of 3, then with x[1] = 2*S/3 the equilibrium
is immediately attained and x[n] = 2*S/3 forever after, since
R is now calculating with integers. E.g. try the above with S<-6
That is what arithmetic ought to be like! But for S not a multiple
of 3 one can get the impression that R is on some sort of drug!

For other values of S (but not all) we observe the same collapse
to x=0 as before, and again it takes 54 steps (ending with x[55]).
Try e.g. S <- 16

For some values of S, however, the iteration ends up in a periodic loop.

For example, with S<-7, at x[52] we get x[52]=4, x[53]=6, x[54]=2,
and then 4 6 2 4 6 2 4 6 2 ... forever (or until Nits cuts in),
so period = 3.

For S<-11, x[52]=8 then 6 then 10 then 2 then 4 then 8 6 10 2 4 ...
so period = 5.

For S<-13, x[51]=4 then 8 10 6 12 2 4 8 10 6 12 2 4 8 ...
so period = 6.

For S<-19, x[51]=12 then 14 10 18 2 4 8 16 6 12 ...
so period = 9.

And so on ...

So, one sniff of something like S<-19, and R is off its head!

All it has to do is multiply by 2 -- and it gets it cumulatively wrong!
R just doesn't add up ...

Season's Greetings to all -- and may your calculations always
be accurate -- to within machine precision ...

Ted.

-------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at wlandres.net>
Date: 22-Dec-2013  Time: 09:59:00
This message was sent by XFMail