[R] bug?

[(Ted Harding)]

On 14-Jul-03 Barry Rowlingson wrote:
>   Of course the main problem is that us 'old-timers' who grew up 
> deciding whether to use REAL*4 or REAL*8 understand what is going on 
> 'under the hood' with floating point numbers, but today's computers are
> being used by mathematicians and statisticians whose perception of 
> numbers comes from a different direction to computer people. "0.7 + 0.1
> does not equal 0.8? Does 1+1 not equal 2? That's madness!", they cry.
> 
>   The solution is education: I'm hoping to teach all our postgraduates
> a short course on computer history and culture, from binary numbers
> up to software engineering, hardware, programming languages etc etc.
> Floating point formats are in there somewhere, for sure.

Spot on, Baz!

An intriguing example to set people is the following:

Write a program to generate a sequence of numbers x(0) ( = anything in
the range 0 <= x <= 1, chosen by user), x(1), x(2), ... , x(n), ...
where
   x(n+1) = f(x(n)), ...
and f() is defined by
   f(x) = 2*x if 0 <= x <= 1/2
   f(x) = 2*(1-x) if 1/2 <= x <= 1

Investigate this series both mathematically and by simulation.

1. Find by mathematics the period-1 equilibrium points [ x(n+1) = x(n) ].

Ans: Clearly: x=0 and x= 2/3

Test: verify these computationally. OOPS!!

2. Find by mathematics the period-2 equilibria [ x(n+2) = x(n) but not
   period-1 ]:

Ans: 2/5 and 4/5

Test: --> OOPS!!

...

k. Find the period-k equilibria

Ans: {2^1, 2^2, ... , 2^k}/(2^k + 1)

Test: --> OOPS!!

And so on. The mathematical equilibria form a set everywhere dense on
[0,1], yet apart from 0 none of them will show up as such on a computer
which is doing standard floating-point arithmetic. (It's a different
story with programs which do exact arithmetic on fractions, preserving
integer numerator and denominator.)

Starting values of the form x(0) = s/2^k head both mathematically and
computationally to 0, arriving in at most (k+1) steps.

Mathematically, all other numbers are not equilibria and the sequence
is "chaotic", yet this cannot be observed either. In fact the sequence
generated will arrive at 0 (and stay there) in at most (N+1) steps where
N is the number of bits in the abscissa of the floating-point
representation.

Every single remotely interesting property of this sequence is inaccesible
by standard arithmetic computation.

Example:
>
f<-function(x){ if((x==0)|(x==1)) stop();
                if(x<=1/2){y<-2*x} else {y<-2*(1-x)}; y }
> x<-2/3; z<-x; while(TRUE){ x<-f(x); z<-c(z,x) }
and then inspect z.

[I invented this example some time ago when bothered by some features of
 a simulated chaos; feel free to use it if it will help to bring home
 some facts of computational life!]

Ted.


--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 167 1972
Date: 14-Jul-03                                       Time: 13:44:45
------------------------------ XFMail ------------------------------