Hi all:

A few days ago, Csaba posted:

*"what was the program, what used by Feigenbaum originally(!), to calculate this number on his HP65? [...] The greatest problem, how to isolate the bifurcation points on the diagram..."*

Actually, it's quite simple. Have a look at this small 4-line, 174-byte *ad-hoc* program I've written for the HP-71B (which can be very easily ported to most any other HP calc, Feigenbaum's own HP-65 included):

10 DESTROY ALL @ DELAY 0,0 @ F=3.2 @ M=0 @ N=1

20 FOR K=2 TO 9 @ P=N+(N-M)/F @ X=FNROOT(P,P,FNB(2^K,FVAR)) @ F=(N-M)/(X-N)

30 PRINT "K=";K;"F=";F @ M=N @ N=X @ NEXT K @ END

40 DEF FNB(K,A) @ B=0 @ FOR I=1 TO K @ B=A-B*B @ NEXT I @ FNB=B @ END DEF

When run, it will perform 9 iterations of a procedure designed to calculate an increasingly refined value of Feigenbaum's constant, starting from a rough initial approximation (3.2). Just press [RUN] and you'll get

these results (K = #iteration, F = estimated value for

Feigenbaum's constant):

K FAfter some 200 seconds, the last result at the end of the 9th iteration is

-------------------------

2 3.21851142201

3 4.38567759894

4 4.60094927695

5 4.65513048088

6 4.66611150124

7 4.66855134603

8 4.66905077160

9 4.66921465563

**F = 4.66921+**, correct to 6 digits within 1 ulp (err = 0.00001+). For comparison purposes, the value of Feigenbaum's constant is

*4.66920 16091029906718+*.

Within the constraints of HP-71B's 12-digit precision and speed, this is more or less the best result you can get. A further 10th iteration would take much, much longer for no real improvement, as the running time grows exponentially with every iteration (there are a number of easy ways of speeding my 4-liner, of course, but I prefer to keep it at its simplest for didactic purposes, if nothing else).

By the time it reaches the 9th iteration, it is finding a root of a *512th-degree* (!) polynomial, which takes a lot of time and incurs in growing cumulative rounding errors when evaluating it a number of times as part of the search for its root. A 10th iteration would mean solving (i.e.: evaluating repeatedly) a 1,024th-degree (!!) polynomial, and cumulative rounding errors mean 6 correct digits is the best you can do in 12-digit arithmetic.

However, if you have access to an arbitrary precision package either for a PC (Maple, Mathematica, etc) or even your HP-71B, you can easily get more correct digits using this algorithm, albeit at a heavy time penalty.

Best regards from V.

*Edited: 3 June 2004, 6:17 a.m. *