[R] bug?
(Ted Harding)
Ted.Harding at nessie.mcc.ac.uk
Mon Jul 14 14:44:45 CEST 2003
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 ------------------------------
More information about the R-help
mailing list