Feigenbaum constant calculation
- Task
Calculate the Feigenbaum constant.
- See
-
- Details in the Wikipedia article: Feigenbaum constant.
AWK
<lang AWK>
- syntax: GAWK -f FEIGENBAUM_CONSTANT_CALCULATION.AWK
BEGIN {
a1 = 1 a2 = 0 d1 = 3.2 max_i = 13 max_j = 10 print(" i d") for (i=2; i<=max_i; i++) { a = a1 + (a1 - a2) / d1 for (j=1; j<=max_j; j++) { x = y = 0 for (k=1; k<=2^i; k++) { y = 1 - 2 * y * x x = a - x * x } a -= x / y } d = (a1 - a2) / (a - a1) printf("%2d %.8f\n",i,d) d1 = d a2 = a1 a1 = a } exit(0)
} </lang>
- Output:
i d 2 3.21851142 3 4.38567760 4 4.60094928 5 4.65513050 6 4.66611195 7 4.66854858 8 4.66906066 9 4.66917155 10 4.66919515 11 4.66920026 12 4.66920098 13 4.66920537
C
<lang c>#include <stdio.h>
void feigenbaum() {
int i, j, k, max_it = 13, max_it_j = 10; double a, x, y, d, a1 = 1.0, a2 = 0.0, d1 = 3.2; printf(" i d\n"); for (i = 2; i <= max_it; ++i) { a = a1 + (a1 - a2) / d1; for (j = 1; j <= max_it_j; ++j) { x = 0.0; y = 0.0; for (k = 1; k <= 1 << i; ++k) { y = 1.0 - 2.0 * y * x; x = a - x * x; } a -= x / y; } d = (a1 - a2) / (a - a1); printf("%2d %.8f\n", i, d); d1 = d; a2 = a1; a1 = a; }
}
int main() {
feigenbaum(); return 0;
}</lang>
- Output:
i d 2 3.21851142 3 4.38567760 4 4.60094928 5 4.65513050 6 4.66611195 7 4.66854858 8 4.66906066 9 4.66917155 10 4.66919515 11 4.66920026 12 4.66920098 13 4.66920537
C++
<lang cpp>#include <iostream>
int main() {
const int max_it = 13; const int max_it_j = 10; double a1 = 1.0, a2 = 0.0, d1 = 3.2;
std::cout << " i d\n"; for (int i = 2; i <= max_it; ++i) { double a = a1 + (a1 - a2) / d1; for (int j = 1; j <= max_it_j; ++j) { double x = 0.0; double y = 0.0; for (int k = 1; k <= 1 << i; ++k) { y = 1.0 - 2.0*y*x; x = a - x * x; } a -= x / y; } double d = (a1 - a2) / (a - a1); printf("%2d %.8f\n", i, d); d1 = d; a2 = a1; a1 = a; }
return 0;
}</lang>
- Output:
i d 2 3.21851142 3 4.38567760 4 4.60094928 5 4.65513050 6 4.66611195 7 4.66854858 8 4.66906066 9 4.66917155 10 4.66919515 11 4.66920026 12 4.66920098 13 4.66920537
C#
<lang csharp>using System;
namespace FeigenbaumConstant {
class Program { static void Main(string[] args) { var maxIt = 13; var maxItJ = 10; var a1 = 1.0; var a2 = 0.0; var d1 = 3.2; Console.WriteLine(" i d"); for (int i = 2; i <= maxIt; i++) { var a = a1 + (a1 - a2) / d1; for (int j = 1; j <= maxItJ; j++) { var x = 0.0; var y = 0.0; for (int k = 1; k <= 1<<i; k++) { y = 1.0 - 2.0 * y * x; x = a - x * x; } a -= x / y; } var d = (a1 - a2) / (a - a1); Console.WriteLine("{0,2:d} {1:f8}", i, d); d1 = d; a2 = a1; a1 = a; } } }
}</lang>
- Output:
i d 2 3.21851142 3 4.38567760 4 4.60094928 5 4.65513050 6 4.66611195 7 4.66854858 8 4.66906066 9 4.66917155 10 4.66919515 11 4.66920026 12 4.66920098 13 4.66920537
D
<lang d>import std.stdio;
void main() {
int max_it = 13; int max_it_j = 10; double a1 = 1.0; double a2 = 0.0; double d1 = 3.2; double a;
writeln(" i d"); for (int i=2; i<=max_it; i++) { a = a1 + (a1 - a2) / d1; for (int j=1; j<=max_it_j; j++) { double x = 0.0; double y = 0.0; for (int k=1; k <= 1<<i; k++) { y = 1.0 - 2.0 * y * x; x = a - x * x; } a -= x / y; } double d = (a1 - a2) / (a - a1); writefln("%2d %.8f", i, d); d1 = d; a2 = a1; a1 = a; }
}</lang>
- Output:
i d 2 3.21851142 3 4.38567760 4 4.60094928 5 4.65513050 6 4.66611195 7 4.66854858 8 4.66906066 9 4.66917155 10 4.66919515 11 4.66920028 12 4.66920099 13 4.66920555
FreeBASIC
<lang freebasic>' version 25-0-2019 ' compile with: fbc -s console
Dim As UInteger i, j, k, maxit = 13, maxitj = 13 Dim As Double x, y, a, a1 = 1, a2, d, d1 = 3.2
Print "Feigenbaum constant calculation:" Print Print " i d" Print "==================="
For i = 2 To maxIt
a = a1 + (a1 - a2) / d1 For j = 1 To maxItJ x = 0 : y = 0 For k = 1 To 2 ^ i y = 1 - 2 * y * x x = a - x * x Next a = a - x / y Next d = (a1 - a2) / (a - a1) Print Using "### ##.#########"; i; d d1 = d a2 = a1 a1 = a
Next
' empty keyboard buffer While Inkey <> "" : Wend Print : Print "hit any key to end program" Sleep End</lang>
- Output:
Feigenbaum constant calculation: i d =================== 2 3.218511422 3 4.385677599 4 4.600949277 5 4.655130495 6 4.666111948 7 4.668548581 8 4.669060660 9 4.669171555 10 4.669195148 11 4.669200285 12 4.669201301 13 4.669198656
Go
<lang go>package main
import "fmt"
func feigenbaum() {
maxIt, maxItJ := 13, 10 a1, a2, d1 := 1.0, 0.0, 3.2 fmt.Println(" i d") for i := 2; i <= maxIt; i++ { a := a1 + (a1-a2)/d1 for j := 1; j <= maxItJ; j++ { x, y := 0.0, 0.0 for k := 1; k <= 1<<uint(i); k++ { y = 1.0 - 2.0*y*x x = a - x*x } a -= x / y } d := (a1 - a2) / (a - a1) fmt.Printf("%2d %.8f\n", i, d) d1, a2, a1 = d, a1, a }
}
func main() {
feigenbaum()
}</lang>
- Output:
i d 2 3.21851142 3 4.38567760 4 4.60094928 5 4.65513050 6 4.66611195 7 4.66854858 8 4.66906066 9 4.66917155 10 4.66919515 11 4.66920026 12 4.66920098 13 4.66920537
Haskell
<lang haskell>import Data.List (mapAccumL)
feigenbaumApprox :: Int -> [Double] feigenbaumApprox mx = snd $ mitch mx 10
where mitch :: Int -> Int -> ((Double, Double, Double), [Double]) mitch mx mxj = mapAccumL (\(a1, a2, d1) i -> let a = iterate (\a -> let (x, y) = iterate (\(x, y) -> (a - (x * x), 1.0 - ((2.0 * x) * y))) (0.0, 0.0) !! (2 ^ i) in a - (x / y)) (a1 + (a1 - a2) / d1) !! mxj d = (a1 - a2) / (a - a1) in ((a, a1, d), d)) (1.0, 0.0, 3.2) [2 .. (1 + mx)]
-- TEST ------------------------------------------------------------------ main :: IO () main =
(putStrLn . unlines) $ zipWith (\i s -> justifyRight 2 ' ' (show i) ++ '\t' : s) [1 ..] (show <$> feigenbaumApprox 13) where justifyRight n c s = drop (length s) (replicate n c ++ s)</lang>
- Output:
1 3.2185114220380866 2 4.3856775985683365 3 4.600949276538056 4 4.6551304953919646 5 4.666111947822846 6 4.668548581451485 7 4.66906066077106 8 4.669171554514976 9 4.669195154039278 10 4.669200256503637 11 4.669200975097843 12 4.669205372040318 13 4.669207514010413
Java
<lang java>public class Feigenbaum {
public static void main(String[] args) { int max_it = 13; int max_it_j = 10; double a1 = 1.0; double a2 = 0.0; double d1 = 3.2; double a;
System.out.println(" i d"); for (int i = 2; i <= max_it; i++) { a = a1 + (a1 - a2) / d1; for (int j = 0; j < max_it_j; j++) { double x = 0.0; double y = 0.0; for (int k = 0; k < 1 << i; k++) { y = 1.0 - 2.0 * y * x; x = a - x * x; } a -= x / y; } double d = (a1 - a2) / (a - a1); System.out.printf("%2d %.8f\n", i, d); d1 = d; a2 = a1; a1 = a; } }
}</lang>
- Output:
i d 2 3.21851142 3 4.38567760 4 4.60094928 5 4.65513050 6 4.66611195 7 4.66854858 8 4.66906066 9 4.66917155 10 4.66919515 11 4.66920026 12 4.66920098 13 4.66920537
Kotlin
<lang scala>// Version 1.2.40
fun feigenbaum() {
val maxIt = 13 val maxItJ = 10 var a1 = 1.0 var a2 = 0.0 var d1 = 3.2 println(" i d") for (i in 2..maxIt) { var a = a1 + (a1 - a2) / d1 for (j in 1..maxItJ) { var x = 0.0 var y = 0.0 for (k in 1..(1 shl i)) { y = 1.0 - 2.0 * y * x x = a - x * x } a -= x / y } val d = (a1 - a2) / (a - a1) println("%2d %.8f".format(i,d)) d1 = d a2 = a1 a1 = a }
}
fun main(args: Array<String>) {
feigenbaum()
}</lang>
- Output:
i d 2 3.21851142 3 4.38567760 4 4.60094928 5 4.65513050 6 4.66611195 7 4.66854858 8 4.66906066 9 4.66917155 10 4.66919515 11 4.66920026 12 4.66920098 13 4.66920537
Modula-2
<lang modula2>MODULE Feigenbaum; FROM FormatString IMPORT FormatString; FROM LongStr IMPORT RealToStr; FROM Terminal IMPORT WriteString,WriteLn,ReadChar;
VAR
buf : ARRAY[0..63] OF CHAR; i,j,k,max_it,max_it_j : INTEGER; a,x,y,d,a1,a2,d1 : LONGREAL;
BEGIN
max_it := 13; max_it_j := 10;
a1 := 1.0; a2 := 0.0; d1 := 3.2;
WriteString(" i d"); WriteLn; FOR i:=2 TO max_it DO a := a1 + (a1 - a2) / d1; FOR j:=1 TO max_it_j DO x := 0.0; y := 0.0; FOR k:=1 TO INT(1 SHL i) DO y := 1.0 - 2.0 * y * x; x := a - x * x END; a := a - x / y END; d := (a1 - a2) / (a - a1); FormatString("%2i ", buf, i); WriteString(buf); RealToStr(d, buf); WriteString(buf); WriteLn; d1 := d; a2 := a1; a1 := a END;
ReadChar
END Feigenbaum.</lang>
Perl
<lang perl>$a1 = 1.0; $a2 = 0.0; $d1 = 3.2;
print " i δ\n";
for my $i (2..13) {
my $a = $a1 + ($a1 - $a2)/$d1; for (1..10) { my $x = 0; my $y = 0; for (1 .. 2**$i) { $y = 1 - 2 * $y * $x; $x = $a - $x*$x; } $a -= $x/$y; }
$d1 = ($a1 - $a2) / ($a - $a1); ($a2, $a1) = ($a1, $a); printf "%2d %17.14f\n", $i, $d1;
}</lang>
- Output:
2 3.21851142203809 3 4.38567759856834 4 4.60094927653806 5 4.65513049539196 6 4.66611194782285 7 4.66854858145148 8 4.66906066077106 9 4.66917155451498 10 4.66919515403928 11 4.66920025650364 12 4.66920097509784 13 4.66920537204032
Perl 6
<lang perl6>my $a1 = 1; my $a2 = 0; my $d = 3.2;
say ' i d';
for 2 .. 13 -> $exp {
my $a = $a1 + ($a1 - $a2) / $d; do { my $x = 0; my $y = 0; for ^2 ** $exp { $y = 1 - 2 * $y * $x; $x = $a - $x²; } $a -= $x / $y; } xx 10; $d = ($a1 - $a2) / ($a - $a1); ($a2, $a1) = ($a1, $a); printf "%2d %.8f\n", $exp, $d;
}</lang>
- Output:
i d 2 3.21851142 3 4.38567760 4 4.60094928 5 4.65513050 6 4.66611195 7 4.66854858 8 4.66906066 9 4.66917155 10 4.66919515 11 4.66920026 12 4.66920098 13 4.66920537
Phix
<lang Phix>constant maxIt = 13,
maxItJ = 10
atom a1 = 1.0,
a2 = 0.0, d1 = 3.2
puts(1," i d\n") for i=2 to maxIt do
atom a = a1 + (a1 - a2) / d1 for j=1 to maxItJ do atom x = 0, y = 0 for k=1 to power(2,i) do y = 1 - 2*y*x x = a - x*x end for a = a - x/y end for atom d = (a1-a2)/(a-a1) printf(1,"%2d %.8f\n",{i,d}) d1 = d a2 = a1 a1 = a
end for</lang>
- Output:
i d 2 3.21851142 3 4.38567760 4 4.60094928 5 4.65513050 6 4.66611195 7 4.66854858 8 4.66906066 9 4.66917155 10 4.66919515 11 4.66920026 12 4.66920098 13 4.66920537
Python
<lang python>max_it = 13 max_it_j = 10 a1 = 1.0 a2 = 0.0 d1 = 3.2 a = 0.0
print " i d" for i in range(2, max_it + 1):
a = a1 + (a1 - a2) / d1 for j in range(1, max_it_j + 1): x = 0.0 y = 0.0 for k in range(1, (1 << i) + 1): y = 1.0 - 2.0 * y * x x = a - x * x a = a - x / y d = (a1 - a2) / (a - a1) print("{0:2d} {1:.8f}".format(i, d)) d1 = d a2 = a1 a1 = a</lang>
- Output:
i d 2 3.21851142 3 4.38567760 4 4.60094928 5 4.65513050 6 4.66611195 7 4.66854858 8 4.66906066 9 4.66917155 10 4.66919515 11 4.66920026 12 4.66920098 13 4.66920537
Racket
<lang racket>#lang racket (define (feigenbaum #:max-it (max-it 13) #:max-it-j (max-it-j 10))
(displayln " i d" (current-error-port)) (define-values (_a _a1 d) (for/fold ((a 1) (a1 0) (d 3.2)) ((i (in-range 2 (add1 max-it)))) (let* ((a′ (for/fold ((a (+ a (/ (- a a1) d)))) ((j (in-range max-it-j))) (let-values (([x y] (for/fold ((x 0) (y 0)) ((k (expt 2 i))) (values (- a (* x x)) (- 1 (* 2 y x)))))) (- a (/ x y))))) (d′ (/ (- a a1) (- a′ a)))) (eprintf "~a ~a\n" (~a i #:width 2) (real->decimal-string d′ 8)) (values a′ a d′)))) d)
(module+ main
(feigenbaum))</lang>
- Output:
i d 2 3.21851142 3 4.38567760 4 4.60094928 5 4.65513050 6 4.66611195 7 4.66854858 8 4.66906066 9 4.66917155 10 4.66919515 11 4.66920026 12 4.66920098 13 4.66920537 4.669205372040318
REXX
<lang rexx>/*REXX pgm calculates the (Mitchell) Feigenbaum bifurcation velocity, #digs can be given*/ parse arg digs maxi maxj . /*obtain optional argument from the CL.*/ if digs== | digs=="," then digs= 30 /*Not specified? Then use the default.*/ if maxi== | maxi=="," then maxi= 20 /* " " " " " " */ if maxJ== | maxJ=="," then maxJ= 10 /* " " " " " " */
- = 4.669201609102990671853203820466201617258185577475768632745651343004134330211314737138,
|| 68974402394801381716 /*◄──Feigenbaum's constant, true value.*/
numeric digits digs /*use the specified # of decimal digits*/
a1= 1 a2= 0 d1= 3.2
say 'Using ' maxJ " iterations for maxJ, with " digs ' decimal digits:' say say copies(' ', 9) center('correct', 11) copies(' ', digs+1) say center('i', 9, "─") center('digits' , 11, '─') center('d', digs+1, "─")
do i=2 for maxi-1 a= a1 + (a1 - a2) / d1 do maxJ x= 0; y= 0 do 2**i; y= 1 - 2 * x * y x= a - x*x end /*2**i*/ a= a - x / y end /*maxj*/ d= (a1 - a2) / (a - a1) /*compute the delta (D) of the function*/ t= max(0, compare(d, #) - 2) /*# true digs so far, ignore dec. point*/ say center(i, 9) center(t, 11) d /*display values for I & D ──►terminal*/ parse value d a1 a with d1 a2 a1 /*assign 3 variables with 3 new values.*/ end /*i*/
say /*stick a fork in it, we're all done. */ say ' true value= ' # / 1 /*true value of Feigenbaum's constant. */</lang>
- output when using the default inputs:
Using 10 iterations for maxJ, with 30 decimal digits: correct ────i──── ──digits─── ───────────────d─────────────── 2 0 3.21851142203808791227050453077 3 1 4.3856775985683390857449485682 4 2 4.60094927653807535781169469969 5 2 4.65513049539198013648625498649 6 3 4.66611194782857138833121364654 7 3 4.66854858144684094804454708811 8 4 4.66906066064826823913257549468 9 4 4.6691715553795113888859465442 10 4 4.66919515603001717402161720542 11 6 4.66920022908685649793393149233 12 7 4.66920131329420417113719511412 13 7 4.66920154578090670783369507315 14 7 4.66920159553749390966169074155 15 9 4.66920160619815215840788706632 16 9 4.66920160848080435144581223484 17 9 4.66920160896974538458267849027 18 10 4.66920160907444981238909862845 19 10 4.66920160909687888294310165196 20 12 4.66920160910169069039564432665 true value= 4.66920160910299067185320382047
Ring
<lang ring>
- Project : Feigenbaum constant calculation
decimals(8) see "Feigenbaum constant calculation:" + nl maxIt = 13 maxItJ = 10 a1 = 1.0 a2 = 0.0 d1 = 3.2 see "i " + "d" + nl for i = 2 to maxIt
a = a1 + (a1 - a2) / d1 for j = 1 to maxItJ x = 0 y = 0 for k = 1 to pow(2,i) y = 1 - 2 * y * x x = a - x * x next a = a - x / y next d = (a1 - a2) / (a - a1) if i < 10 see "" + i + " " + d + nl else see "" + i + " " + d + nl ok d1 = d a2 = a1 a1 = a
next </lang> Output:
Feigenbaum constant calculation: i d 2 3.21851142 3 4.38567760 4 4.60094928 5 4.65513050 6 4.66611195 7 4.66854858 8 4.66906066 9 4.66917155 10 4.66919515 11 4.66920026 12 4.66920098 13 4.66920537
Sidef
<lang ruby>var a1 = 1 var a2 = 0 var δ = 3.2.float
say " i\tδ"
for i in (2..15) {
var a0 = ((a1 - a2)/δ + a1) 10.times { var (x, y) = (0, 0) 2**i -> times { y = (1 - 2*x*y) x = (a0 - x²) } a0 -= x/y } δ = ((a1 - a2) / (a0 - a1)) (a2, a1) = (a1, a0) printf("%2d %.8f\n", i, δ)
}</lang>
- Output:
i δ 2 3.21851142 3 4.38567760 4 4.60094928 5 4.65513050 6 4.66611195 7 4.66854858 8 4.66906066 9 4.66917156 10 4.66919516 11 4.66920023 12 4.66920131 13 4.66920155 14 4.66920160 15 4.66920161
zkl
<lang zkl>fcn feigenbaum{
maxIt,maxItJ,a1,a2,d1,a,d := 13, 10, 1.0, 0.0, 3.2, 0, 0; println(" i d"); foreach i in ([2..maxIt]){ a=a1 + (a1 - a2)/d1; foreach j in ([1..maxItJ]){ x,y := 0.0, 0.0;
foreach k in ([1..(1).shiftLeft(i)]){ y,x = 1.0 - 2.0*y*x, a - x*x; } a-=x/y
} d=(a1 - a2)/(a - a1); println("%2d %.8f".fmt(i,d)); d1,a2,a1 = d,a1,a; }
}();</lang>
- Output:
i d 2 3.21851142 3 4.38567760 4 4.60094928 5 4.65513050 6 4.66611195 7 4.66854858 8 4.66906066 9 4.66917155 10 4.66919515 11 4.66920026 12 4.66920098 13 4.66920537