Feigenbaum constant calculation

From Rosetta Code
Feigenbaum constant calculation is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.


Task

Calculate the Feigenbaum constant.


See



AWK[edit]

 
# 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)
}
 
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[edit]

Translation of: Ring
#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;
}
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#[edit]

Translation of: Kotlin
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;
}
}
}
}
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[edit]

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;
}
}
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

Go[edit]

Translation of: Ring
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()
}
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[edit]

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)
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[edit]

Translation of: Kotlin
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;
}
}
}
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[edit]

Translation of: Ring
// 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()
}
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[edit]

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.

Perl[edit]

$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;
}
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[edit]

Works with: Rakudo version 2018.04.01
Translation of: Ring
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;
}
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[edit]

Translation of: Ring
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
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[edit]

Translation of: C
#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))
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[edit]

Translation of: Sidef
/*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. */
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[edit]

 
# 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
 

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[edit]

Translation of: Perl 6
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 -)
}
a0 -= x/y
}
δ = ((a1 - a2) / (a0 - a1))
(a2, a1) = (a1, a0)
printf("%2d %.8f\n", i, δ)
}
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[edit]

Translation of: Kotlin
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;
}
}();
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