Re: [R] Season's Greetings (and great news ... )!
"(or whatever the keyboard analogue of that may be) "
Hands clasped? Fingers interlaced?
John Kane
Kingston ON Canada
> -Original Message-
> From: ted.hard...@wlandres.net
> Sent: Sun, 22 Dec 2013 18:37:18 - (GMT)
> To: r-help@r-project.org
> Subject: Re: [R] Season's Greetings (and great news ... )!
>
> Thanks for the comments, Bert and Mehmet! It is of course a serious
> and interesting area to explore (and I'm aware of the chaos context;
> I initially got into this areas year ago when I was exploring the
> possibilities for chaos in fish population dynamics -- and they're
> certainly there)!
>
> But, before anyone takes my posting *too* seriously, let me say that
> it was written tongue-in-cheek (or whatever the keyboard analogue of
> that may be). I'm certainly not "blaming R".
>
> Have fun anyway!
> Ted.
>
> On 22-Dec-2013 17:35:56 Bert Gunter wrote:
>> Yes.
>>
>> See also Feigenbaum's constant and chaos theory for the general context.
>>
>> Cheers,
>> Bert
>>
>> On Sun, Dec 22, 2013 at 8:54 AM, Suzen, Mehmet wrote:
>>> 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
>>> 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.7
>>>> 2: 0.7
>>>> 3: 0.7
>>>> 4: 0.7
>>>> 5: 0.7
>>>> ...
>>>>
>>>> but, later on,
>>>>
>>>> 24: 0.7
>>>> 25: 0.6
>>>> 26: 0.8
>>>> 2
Re: [R] Season's Greetings (and great news ... )!
Thanks for the comments, Bert and Mehmet! It is of course a serious
and interesting area to explore (and I'm aware of the chaos context;
I initially got into this areas year ago when I was exploring the
possibilities for chaos in fish population dynamics -- and they're
certainly there)!
But, before anyone takes my posting *too* seriously, let me say that
it was written tongue-in-cheek (or whatever the keyboard analogue of
that may be). I'm certainly not "blaming R".
Have fun anyway!
Ted.
On 22-Dec-2013 17:35:56 Bert Gunter wrote:
> Yes.
>
> See also Feigenbaum's constant and chaos theory for the general context.
>
> Cheers,
> Bert
>
> On Sun, Dec 22, 2013 at 8:54 AM, Suzen, Mehmet wrote:
>> 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 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.7
>>> 2: 0.7
>>> 3: 0.7
>>> 4: 0.7
>>> 5: 0.7
>>> ...
>>>
>>> but, later on,
>>>
>>> 24: 0.7
>>> 25: 0.6
>>> 26: 0.8
>>> 27: 0.4
>>> 28: 0.66672
>>> ...
>>>
>>> 46: 0.667968750
>>> 47: 0.664062500
>>> 48: 0.671875000
>>> 49: 0.65625
>>> 50: 0.68750
>>> 51: 0.62500
>>> 52: 0.75000
>>> 53: 0.5
>>> 54: 1.0
>>> 55: 0.0
>>>
>>> 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.000 -> 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.101 x[5] -> 2*(1 - 0.101) = 2*0.011 ->
>>> 0.110 x[6] -> 2*(1 - 0.110) = 2*0.010 ->
>>> 0.100 x[7] -> 2*0.100 = 1.000 ->
>>> 1.000 x[8] -> 2*(1 - 1.000) = 2*0 ->
>>> 0.000 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
>>>
Re: [R] Season's Greetings (and great news ... )!
Yes.
See also Feigenbaum's constant and chaos theory for the general context.
Cheers,
Bert
On Sun, Dec 22, 2013 at 8:54 AM, Suzen, Mehmet wrote:
> 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 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.7
>> 2: 0.7
>> 3: 0.7
>> 4: 0.7
>> 5: 0.7
>> ...
>>
>> but, later on,
>>
>> 24: 0.7
>> 25: 0.6
>> 26: 0.8
>> 27: 0.4
>> 28: 0.66672
>> ...
>>
>> 46: 0.667968750
>> 47: 0.664062500
>> 48: 0.671875000
>> 49: 0.65625
>> 50: 0.68750
>> 51: 0.62500
>> 52: 0.75000
>> 53: 0.5
>> 54: 1.0
>> 55: 0.0
>>
>> 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.000 -> 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.101 x[5] -> 2*(1 - 0.101) = 2*0.011 ->
>> 0.110 x[6] -> 2*(1 - 0.110) = 2*0.010 ->
>> 0.100 x[7] -> 2*0.100 = 1.000 ->
>> 1.000 x[8] -> 2*(1 - 1.000) = 2*0 ->
>> 0.000 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
Re: [R] 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 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.7
> 2: 0.7
> 3: 0.7
> 4: 0.7
> 5: 0.7
> ...
>
> but, later on,
>
> 24: 0.7
> 25: 0.6
> 26: 0.8
> 27: 0.4
> 28: 0.66672
> ...
>
> 46: 0.667968750
> 47: 0.664062500
> 48: 0.671875000
> 49: 0.65625
> 50: 0.68750
> 51: 0.62500
> 52: 0.75000
> 53: 0.5
> 54: 1.0
> 55: 0.0
>
> 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.000 -> 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.101 x[5] -> 2*(1 - 0.101) = 2*0.011 ->
> 0.110 x[6] -> 2*(1 - 0.110) = 2*0.010 ->
> 0.100 x[7] -> 2*0.100 = 1.000 ->
> 1.000 x[8] -> 2*(1 - 1.000) = 2*0 ->
> 0.000 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 va
[R] Season's Greetings (and great news ... )!
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.7
2: 0.7
3: 0.7
4: 0.7
5: 0.7
...
but, later on,
24: 0.7
25: 0.6
26: 0.8
27: 0.4
28: 0.66672
...
46: 0.667968750
47: 0.664062500
48: 0.671875000
49: 0.65625
50: 0.68750
51: 0.62500
52: 0.75000
53: 0.5
54: 1.0
55: 0.0
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.000 -> 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.101 x[5] -> 2*(1 - 0.101) = 2*0.011 ->
0.110 x[6] -> 2*(1 - 0.110) = 2*0.010 ->
0.100 x[7] -> 2*0.100 = 1.000 ->
1.000 x[8] -> 2*(1 - 1.000) = 2*0 ->
0.000 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
