🚀 go-pugleaf

/*
LYNCH%XAVIER.BITNET@uga.cc.uga.edu wrote on  Wed, 13 Oct 1993 18:48:00 EST

|> The Feigenbaum number is hailed in the New Guiness book as the
|>         "Newest Mathematical Constant"
|> 
|> Could someone please explain to me what the Feigenbaum number represents?
|> Also, these people producing the number to many digits, do you have
|> at least the pseudo code for this feat? I would appreciate seeing it.


From fractal-faq Oct 2, 1993 by Ken Shirriff

Subject: Feigenbaum's constant

Q10: What is Feigenbaum's constant?
A10: In a period doubling cascade, such as the logistic equation, consider the
parameter values where period-doubling events occur.  The limit of the ratio
of distances between consecutive doubling values is Feigenbaum's constant.
Based on independent computations by Jay Hill and Keith Briggs, it has the
value 4.669201609102990671853...  Note: several books have published incorrect
values starting 4.66920166...; the last repeated 6 is a typographical error.

Here is a C++ program which illustrates the Feigenbaum constant calculation
using the Mandelbrot Set iteration formula.  Most texts use the logistic
formula  x=r*x*(1-x).   This program shows the constant computed for two
period doubling cascades.  One starts with the period 1 cardioid and the
second starts with the period 3 cardioid.  At the end of the cascades are the
Feigenbaum points, for which I compute -1.40115518909205060052382678789... 
and -1.7798180758349913642620671...   Change the comment format and this 
program becomes C.  I prefer C++, a simple change of double to complex
here and there and you can explore the cascades from complex components.
When you do, you can try n=8, b=(-0.156520166833755062,-1.0322471089228318),
a=(-0.015608533223847222 -1.03664566584786215).
*/
// C++ Compute Feigenbaum Universal Constant
// (C) Jay R Hill, 1993

#include <stdlib.h>
double Newton(long n, double c){
// find root of z=z^2 + c iterated n times
  double z, z1, nc=c, absz=1; long i,j;
  j=0;  // j limits iterations
  while(j<  7 && absz>(double)1e-13){ ++j;
    z=(double)0;
    z1=(double)0;
    for(i=0; i<n; ++i){
      z1=2*z1*z+1;        // first derivative
      z=z*z+nc;           // iterate
    }
    nc-=z/z1;         // Newton update
    absz=fabs(z);     // for early out
  }
  return(nc);
}
void go(long n0, double a, double b){
  double F=4.0, Tmp=a,NewC=a+(a-b)/F,OldC=b; long i,n=2*n0;
  for(i=0; i<10; ++i){
    NewC=Newton(n, NewC); 
// make new estimate for Feigenbaum Universal Constant
    F=(Tmp-OldC)/(NewC-Tmp);  
    printf("%.16lf %.16lf %.16lf %.16lf\n",OldC,Tmp,NewC,F);
    OldC=Tmp;
    Tmp=NewC;
// use estimate for Feigenbaum Universal Constant to find next 
// period doubling cascade
    NewC+=(NewC-OldC)/(F);
    n*=2;
  }
}
void main(void)
{
  double a,b; long n;
  printf(
"      c1                   c2                c3             F: (c2-c1)/(c3-c2)\n");
  b=-0.0;
  a=-1.0;
  n=2; 
  go(n,a,b);  // look at sequence starting at period 1 cardioid
  printf("\n");
  b=-1.75487766624669276;
  a=-1.7728929033816238;
  n=6;
  go(n,a,b);  // look at sequence starting at period 3 cardioid
}
/*
      c1                   c2                c3             F: (c2-c1)/(c3-c2)
-0.0000000000000000 -1.0000000000000000 -1.3107026413368328 3.2185114220380888
-1.0000000000000000 -1.3107026413368328 -1.3815474844320612 4.3856775985683472
-1.3107026413368328 -1.3815474844320612 -1.3969453597045607 4.6009492765379756
-1.3815474844320612 -1.3969453597045607 -1.4002530812147831 4.6551304953917194
-1.3969453597045607 -1.4002530812147831 -1.4009619629448422 4.6661119478231585
-1.4002530812147831 -1.4009619629448422 -1.4011138049397749 4.6685485815265810
-1.4009619629448422 -1.4011138049397749 -1.4011463258269472 4.6690606602496461
-1.4011138049397749 -1.4011463258269472 -1.4011532908499256 4.6691715552590010
-1.4011463258269472 -1.4011532908499256 -1.4011547825466131 4.6691951766774649
-1.4011532908499256 -1.4011547825466131 -1.4011551020224637 4.6692001427014445

-1.7548776662466927 -1.7728929033816239 -1.7782668211110817 3.3523470290916921
-1.7728929033816239 -1.7782668211110817 -1.7794823875518131 4.4209164957074565
-1.7782668211110817 -1.7794823875518131 -1.7797460275941450 4.6107049216787335
-1.7794823875518131 -1.7797460275941450 -1.7798026381018630 4.6570866957290749
-1.7797460275941450 -1.7798026381018630 -1.7798147692180581 4.6665539104053879
-1.7798026381018630 -1.7798147692180581 -1.7798173676438036 4.6686407015227296
-1.7798147692180581 -1.7798173676438036 -1.7798179241614378 4.6690807007454609
-1.7798173676438036 -1.7798179241614378 -1.7798180433511361 4.6691756249117882
-1.7798179241614378 -1.7798180433511361 -1.7798180688779419 4.6691975235440442
-1.7798180433511361 -1.7798180688779419 -1.7798180743450029 4.6692008116173804

*/
Warmly,
Jay "Not tonight honey, I'm in a period doubling cascade' Hill
--
{ hilljr@jupiter.saic.com } begin writeln(3*ln(640320)/sqrt(163):17:15) end.
void main(){double sqrt(), y=1/sqrt(2.), a=.5, m=1,z; int n=0;
for(;m*=2,z=sqrt(1-y*y),y=(1-z)/(1+z),a=a*(1+y)*(1+y)-m*y,n<4;n++);
printf("17.15lf\n",1/a);}