Feigenbaum constant calculation

From Rosetta Code
Task
Feigenbaum constant calculation
You are encouraged to solve this task according to the task description, using any language you may know.
Task

Calculate the Feigenbaum constant.


See



11l

Translation of: Python
V max_it = 13
V max_it_j = 10
V a1 = 1.0
V a2 = 0.0
V d1 = 3.2
V a = 0.0

print(‘ i       d’)
L(i) 2..max_it
   a = a1 + (a1 - a2) / d1
   L(j) 1..max_it_j
      V x = 0.0
      V y = 0.0
      L(k) 1..(1 << i)
         y = 1.0 - 2.0 * y * x
         x = a - x * x
      a = a - x / y
   V d = (a1 - a2) / (a - a1)
   print(‘#2    #.8’.format(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

Ada

Translation of: Ring
with Ada.Text_IO;         use Ada.Text_IO;
with Ada.Integer_Text_IO; use Ada.Integer_Text_IO;

procedure Main is
   procedure feigenbaum is
      subtype i_range is Integer range 2 .. 13;
      subtype j_range is Integer range 1 .. 10;

      -- the number of digits in type Real is reduced to 15 to produce the
      -- results reported by C, C++, C# and Ring. Increasing the number of
      -- digits in type Real produces the results reported by D.

      type Real is digits 15;
      package Real_Io is new Float_IO (Real);
      use Real_Io;

      a, x, y, d : Real;
      a1         : Real := 1.0;
      a2         : Real := 0.0;
      d1         : Real := 3.2;
   begin
      Put_Line (" i       d");
      for i in i_range loop
         a := a1 + (a1 - a2) / d1;
         for j in j_range loop
            x := 0.0;
            y := 0.0;
            for k in 1 .. 2**i loop
               y := 1.0 - 2.0 * x * y;
               x := a - x * x;
            end loop;
            a := a - x / y;
         end loop;
         d := (a1 - a2) / (a - a1);
         Put (Item => i, Width => 2);
         Put (Item => d, Fore => 5, Aft => 8, Exp => 0);
         New_Line;
         d1 := d;
         a2 := a1;
         a1 := a;
      end loop;
   end feigenbaum;

begin
   feigenbaum;
end Main;
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

ALGOL 68

Works with: ALGOL 68G version Any - tested with release 2.8.3.win32
Translation of: Ring
# Calculate the Feigenbaum constant #
 
print( ( "Feigenbaum constant calculation:", newline ) );
INT max it   = 13;
INT max it j = 10;
REAL a1 := 1.0;
REAL a2 := 0.0;
REAL d1 := 3.2;
print( ( "i  ", "d", newline ) );
FOR i FROM 2 TO max it DO
     REAL a := a1 + (a1 - a2) / d1;
     FOR j TO max it j DO
          REAL x := 0;
          REAL y := 0;
          FOR k TO 2 ^ i DO
               y := 1 - 2 * y * x;
               x := a - x * x 
          OD;
          a := a - x / y
     OD;
     REAL d = (a1 - a2) / (a - a1);
     IF i < 10 THEN
        print( ( whole( i, 0 ), "  ", fixed( d, -10, 8 ), newline ) )
     ELSE
        print( ( whole( i, 0 ), " ",  fixed( d, -10, 8 ), newline ) )
     FI;
     d1 := d;
     a2 := a1;
     a1 := a
OD
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

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

BASIC

BASIC256

maxIt = 13 : maxItj = 13
a1 = 1.0 : a2 = 0.0 : d = 0.0 : 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.0 : y = 0.0
        for k = 1 to 2 ^ i
            y = 1 - 2 * y * x
            x = a - x * x
        next k
        a -= x / y
    next j
    d = (a1 - a2) / (a - a1)
    print rjust(i,3); chr(9); ljust(d,13,"0")
    d1 = d
    a2 = a1
    a1 = a
next i
Output:
Same as FreeBASIC entry.

Chipmunk Basic

Works with: Chipmunk Basic version 3.6.4
Translation of: FreeBASIC
100 cls
110 mit = 13
120 mitj = 13
130 a1 = 1
140 a2 = 0
150 d = 0
160 d1 = 3.2
170 print "Feigenbaum constant calculation:"
180 print
190 print "  i     d"
200 print "==================="
210 for i = 2 to mit
220   a = a1+(a1-a2)/d1
230   for j = 1 to mitj
240     x = 0
250     y = 0
260     for k = 1 to 2^i
270       y = 1-2*y*x
280       x = a-x*x
290     next k
300     a = a-(x/y)
310   next j
320   d = (a1-a2)/(a-a1)
330   print using "###";i;"    ";
335   print using "##.#########";d
340   d1 = d
350   a2 = a1
360   a1 = a
370 next i
380 end

Just BASIC

maxit = 13 : maxitj = 13
a1 = 1.0 : a2 = 0.0 : d = 0.0 : 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 k
    a = a - (x / y)
  next j
  d = (a1 - a2) / (a - a1)
  print i; tab(8); d
  d1 = d
  a2 = a1
  a1 = a
next i

MSX Basic

Works with: MSX BASIC version any
Translation of: FreeBASIC
100 CLS
110 mit = 13
120 mitj = 13
130 a1 = 1
140 a2 = 0
150 d = 0
160 d1 = 3.2
170 PRINT "Feigenbaum constant calculation:"
180 PRINT
190 PRINT "  i     d"
200 PRINT "==================="
210 FOR i = 2 TO mit
220   a = a1 + (a1 - a2) / d1
230   FOR j = 1 TO mitj
240     x = 0
250     y = 0
260     FOR k = 1 TO 2 ^ i
270       y = 1 - 2 * y * x
280       x = a - x * x
290     NEXT k
300     a = a - (x / y)
310   NEXT j
320   d = (a1 - a2) / (a - a1)
330   PRINT USING "###    ##.#########"; i; d
340   d1 = d
350   a2 = a1
360   a1 = a
370 NEXT i
380 END

True BASIC

LET maxit = 13
LET maxitj = 13
LET a1 = 1.0
LET d1 = 3.2

PRINT "Feigenbaum constant calculation:"
PRINT
PRINT "  i     d"
PRINT "==================="

FOR i = 2 to maxit
    LET a = a1 + (a1 - a2) / d1
    FOR j = 1 to maxitj
        LET x = 0
        LET y = 0
        FOR k = 1 to 2 ^ i
            LET y = 1 - 2 * y * x
            LET x = a - x * x
        NEXT k
        LET a = a - (x / y)
    NEXT j
    LET d = (a1 - a2) / (a - a1)
    PRINT using "###    ##.#########": i, d
    LET d1 = d
    LET a2 = a1
    LET a1= a
NEXT i
END
Output:
Same as FreeBASIC entry.

Yabasic

maxIt = 13 : maxItj = 13
a1 = 1.0 : a2 = 0.0 : d = 0.0 : d1 = 3.2

print "Feigenbaum constant calculation:"
print "\n  i      d"
print "===================="

for i = 2 to maxIt
    a = a1 + (a1 - a2) / d1
    for j = 1 to maxItj
        x = 0.0 : y = 0.0
        for k = 1 to 2 ^ i
            y = 1 - 2 * y * x
            x = a - x * x
        next k
        a = a - x / y
    next j
    d = (a1 - a2) / (a - a1)
    print i using("###"), chr$(9), d
    d1 = d
    a2 = a1
    a1 = a
next i

C

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#

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

C++

Translation of: C
#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;
}
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

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

Delphi

Works with: Delphi version 6.0

Translated from Algol

procedure FeigenbaumConstant(Memo: TMemo);
{ Calculate the Feigenbaum constant }
const IMax = 13;
const JMax = 10;
var I,J,K: integer;
var A1,A2,D1,X,Y: double;
var A,D: double;
begin
Memo.Lines.Add('Feigenbaum constant calculation:');
{Set initial starting values for iterations}
A1:=1.0; A2:=0.0; D1:=3.2;
Memo.Lines.Add(' I           A           D');
for I:=2 to IMax do
	begin
	{Find next Bifurcation parameter, A}
	A:=A1 + (A1 - A2) / D1;
	for J:=1 to JMax do
		begin
		X:=0; Y:=0;
		for K:=1 to 1 shl i do
			begin
			Y:=1 - 2 * y * x;
			X:=A - X * X
			end;
		A:=A - X / Y
		end;
	{Use current and previous values of A}
	{to calculate the Feigenbaum constant D }
	D:= (A1 - A2) / (A - A1);
	Memo.Lines.Add(Format('%2d  %2.8f  %2.8f',[I,A,D]));
	D1:=D; A2:=A1; A1:=A;
	end;
end;
Output:
Feigenbaum constant calculation:
 I           A           D
 2  1.31070264  3.21851142
 3  1.38154748  4.38567760
 4  1.39694536  4.60094928
 5  1.40025308  4.65513050
 6  1.40096196  4.66611195
 7  1.40111380  4.66854858
 8  1.40114633  4.66906066
 9  1.40115329  4.66917155
10  1.40115478  4.66919515
11  1.40115510  4.66920028
12  1.40115517  4.66920099
13  1.40115519  4.66920555

EasyLang

Translation of: AWK
numfmt 6 0
a1 = 1 ; a2 = 0 ; d1 = 3.2
ipow2 = 4
for i = 2 to 13
   a = a1 + (a1 - a2) / d1
   for j = 1 to 10
      x = 0 ; y = 0
      for k = 1 to ipow2
         y = 1 - 2 * y * x
         x = a - x * x
      .
      a -= x / y
   .
   d = (a1 - a2) / (a - a1)
   print i & " " & d
   d1 = d ; a2 = a1 ; a1 = a
   ipow2 *= 2
.

F#

Translation of: C#
open System

[<EntryPoint>]
let main _ = 
    let maxIt = 13
    let maxItJ = 10
    let mutable a1 = 1.0
    let mutable a2 = 0.0
    let mutable d1 = 3.2
    Console.WriteLine(" i       d")
    for i in 2 .. maxIt do
        let mutable a = a1 + (a1 - a2) / d1
        for j in 1 .. maxItJ do
            let mutable x = 0.0
            let mutable y = 0.0
            for _ in 1 .. (1 <<< i) do
                y <- 1.0 - 2.0 * y * x
                x <- a - x * x
            a <- a - x / y
        let d = (a1 - a2) / (a - a1)
        Console.WriteLine("{0,2:d}    {1:f8}", i, d)
        d1 <- d
        a2 <- a1
        a1 <- a
    0 // return an integer exit code
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

Factor

Translation of: Raku
USING: formatting io locals math math.ranges sequences ;

[let
      1 :> a1!
      0 :> a2!
    3.2 :> d!

    " i d" print

    2 13 [a,b] [| exp |
        a1 a2 - d /f a1 + :> a!
        10 [
            0 :> x!
            0 :> y!
            exp 2^ [
                1 2 x y * * - y!
                a x sq - x!
            ] times
            a x y /f - a!
        ] times
        a1 a2 - a a1 - /f d!
        a1 a2! a a1!
        exp d "%2d %.8f\n" printf
    ] each
]
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

Fortran

      program feigenbaum
      implicit none

      integer i, j, k
      real ( KIND = 16 ) x, y, a, b, a1, a2, d1

      print '(a4,a13)', 'i', 'd'

      a1 = 1.0;
      a2 = 0.0;
      d1 = 3.2;

      do i=2,20
         a = a1 + (a1 - a2) / d1;
         do j=1,10
            x = 0
            y = 0
            do k=1,2**i
                y = 1 - 2 * y * x;
                x = a - x**2;
            end do
            a = a - x / y;
         end do

         d1 = (a1 - a2) / (a - a1);
         a2 = a1;
         a1 = a;
         print '(i4,f13.10)', i, d1
     end do
     end
Output:
   i            d
   2 3.2185114220
   3 4.3856775986
   4 4.6009492765
   5 4.6551304954
   6 4.6661119478
   7 4.6685485814
   8 4.6690606606
   9 4.6691715554
  10 4.6691951560
  11 4.6692002291
  12 4.6692013133
  13 4.6692015458
  14 4.6692015955
  15 4.6692016062
  16 4.6692016085
  17 4.6692016090
  18 4.6692016091
  19 4.6692016091
  20 4.6692016091

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

Fōrmulæ

Fōrmulæ programs are not textual, visualization/edition of programs is done showing/manipulating structures but not text. Moreover, there can be multiple visual representations of the same program. Even though it is possible to have textual representation —i.e. XML, JSON— they are intended for storage and transfer purposes more than visualization and edition.

Programs in Fōrmulæ are created/edited online in its website.

In this page you can see and run the program(s) related to this task and their results. You can also change either the programs or the parameters they are called with, for experimentation, but remember that these programs were created with the main purpose of showing a clear solution of the task, and they generally lack any kind of validation.

Solution

Test case

FutureBasic

Translation of: Ring and Phix
window 1, @"Feignenbaum Constant", ( 0, 0, 200, 300 )

_maxIt  = 13
_maxItJ = 10

void local fn Feignenbaum
  NSUInteger i, j, k
  double     a1 = 1.0, a2 = 0.0, d1 = 3.2
  
  print "Feignenbaum Constant"
  print " i       d"
  
  for i = 2 to _maxIt
    double a = a1 + ( a1 - a2 ) / d1
    for j = 1 to _maxItJ
      double x = 0, y = 0
      for k = 1 to fn pow( 2, i )
        y = 1 - 2 * y * x
        x = a - x * x
      next
      a = a - x / y
    next
    double d = ( a1 - a2 ) / ( a - a1 )
    printf @"%2d.   %.8f", i, d
    d1 = d
    a2 = a1
    a1 = a
  next
end fn

fn Feignenbaum

HandleEvents
Output:
 Feignenbaum Constant
 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

Go

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


Groovy

Translation of: Java
class Feigenbaum {
    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

        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)
            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

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)
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

J

Translated from the beautiful Fōrmulæ version. Rather than a verb, the conjunction pre-assigns m and n .

Feigenbaum =: conjunction define  NB. use:  n Feigenbaum m
 irange=: <. + i.@:>:@:|@:-  NB. inclusive range
 a=. 0 1
 delta=. , 3.2
 for_i. 3 irange n do.
  tmp=. ({: + ({:delta) *inv ({: - _2&{)) a
  for. i. m do.
   'b bp'=. 0
   for. i. 2 ^ <: i do.
    'b bp'=. (tmp - *: b) , 1 _2 p. b * bp
   end.
   tmp=. tmp - b % bp
  end.
  a=. a , tmp
  delta=. delta , %/@:(-/"1) (- 2 3 ,: 1 2) { a
 end.
 2 14j6 14j6 ": (#\ i. # delta) ,. (}. a) ,. delta
)

   8  Feigenbaum  13
 1      1.000000      3.200000
 2      1.310703      3.218511
 3      1.381547      4.385678
 4      1.396945      4.600949
 5      1.400253      4.655130
 6      1.400962      4.666112
 7      1.401114      4.668549
 8      1.401146      4.669061
 9      1.401153      4.669172
10      1.401155      4.669195
11      1.401155      4.669200
12      1.401155      4.669201

Java

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

jq

def feigenbaum_delta(imax; jmax):
  def lpad: tostring | (" " *  (4 - length)) + .;
  def pp(i;x): "\(i|lpad)   \(x)";

  "Feigenbaum's delta constant incremental calculation:",
     pp("i"; "δ"),
     pp(1; "3.20"),
     ( foreach range(2; 1+imax) as $i (
         {a1: 1.0, a2: 0.0, d1: 3.2};

         .a = .a1 + (.a1 - .a2) / .d1
         | reduce range(1; 1+jmax) as $j (.;
             .x = 0 | .y = 0
             | reduce range(1; 1+pow(2;$i)) as $k (.;
                 .y = (1 - 2 * .x * .y)
                 | .x = .a - (.x * .x) )
             | .a -= (.x / .y) )
         | .d = (.a1 - .a2) / (.a - .a1)
         | .d1 = .d | .a2 = .a1 | .a1 = .a; 
	    pp($i; .d) ) ) ;
 
Feigenbaum_delta(13; 10)
Output:
Feigenbaum's delta constant incremental calculation:
   i   δ
   1   3.20
   2   3.2185114220380866
   3   4.3856775985683365
   4   4.600949276538056
   5   4.6551304953919646
   6   4.666111947822846
   7   4.668548581451485
   8   4.66906066077106
   9   4.669171554514976
  10   4.669195154039278
  11   4.669200256503637
  12   4.669200975097843
  13   4.669205372040318

Julia

# http://en.wikipedia.org/wiki/Feigenbaum_constant

function feigenbaum_delta(imax=23, jmax=20)
    a1, a2, d1 = BigFloat(1.0), BigFloat(0.0), BigFloat(3.2)
    println("Feigenbaum's delta constant incremental calculation:\ni   δ\n1   3.20")
    for i in 2:imax
        a = a1 + (a1 - a2) / d1
        for j in 1:jmax
            x, y = 0, 0
            for k in 1:2^i
                y = 1 - 2 * x * y
                x = a - x * x
            end
            a -= x / y
        end
        d = (a1 - a2) / (a - a1)
        println(rpad(i, 4), lpad(d, 4))
        d1, a2 = d, a1
        a1 = a
    end
end

feigenbaum_delta()
Output:
Feigenbaum's delta constant incremental calculation:
i   δ
1   3.20
2   3.218511422038087912270504530742813256028820377971082199141994437483271226037533
3   4.385677598568339085744948568775522346103216356576497808699630752612705940390646
4   4.600949276538075357811694698623834985023552496633543372295593454454329771521727
5   4.655130495391980136486254995856898819475460497385226078363311588165123307017281
6   4.66611194782857138833121369671177648071905897173694216397236891198998639455025
7   4.668548581446840948044543680148146265543287896654348757317309551400403337843036
8   4.66906066064826823913259982263027263779968209542149740052288679867743088942764
9   4.669171555379511388886004609897567088240676573170789783804375113804695091803033
10  4.669195156030017174021108801191492093392147908605756405516325961597435372704323
11  4.669200229086856497938353781004067217408888048906823830162962242800074595934665
12  4.669201313294204171164754941185571183728248888986548913352217226469150028661929
13  4.669201545780906707506058109930429736431564330452605295006142805341042630340361
14  4.669201595537493910292470639289646040074547412490596040512777985387237785978782
15  4.669201606198152157723831097078594524421336516011873717994000712976201143278191
16  4.669201608480804423294067945898622842792868381815074127672747764898152898198069
17  4.669201608969744700482485321938373343907385540992447405883605282416375303280911
18  4.669201609074452566227981520370886753946099646679618270214759101315481224820708
19  4.669201609096878794705135037864783677622666525741836726064298799595215295927305
20  4.66920160910168168118696016084580172992808889324407617097679098039831535247408
21  4.669201609102710327837210208629111857781724142614997392167298168695631199065625
22  4.669201609102930630539778141205517641783439121041016813735799961205502985593042
23  4.66920160910297781286849594159066394676896043144121209732784416240857379387701

Kotlin

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

Lambdatalk

Following the Python code in a recursive mode.

{feigenbaum 11}  // on my computer stackoverflow for values greater than 11
-> [3.2185114220380866,4.3856775985683365,4.600949276538056,4.6551304953919646,4.666111947822846,
4.668548581451485,4.66906066077106,4.669171554514976,4.669195154039278,4.669200256503637]

with:

{def feigenbaum
 {lambda {:maxi}
  {f3 :maxi 10 1 0 3.2 0 {A.new} 2}}}

{def f3
 {lambda {:maxi :maxj :a1 :a2 :d1 :a3 :s :i}
  {if {< :i {+ :maxi 1}}
   then {let { {:maxi :maxi} {:maxj :maxj} {:a1 :a1} {:a2 :a2}
               {:a3 {f2 {+ :a1 {/ {- :a1 :a2} :d1}} :i :maxj 1} }
               {:s :s} {:i :i}
             } {f3 :maxi :maxj :a3 :a1 {/ {- :a1 :a2} {- :a3 :a1}} :a3
                   {A.addlast! {/ {- :a1 :a2} {- :a3 :a1}} :s} {+ :i 1}} }
   else :s}}}

{def f2
 {lambda {:a :i :maxj :j}
  {if {< :j {+ :maxj 1}}
   then {f2 {f1 :a :i 0 0 1} :i :maxj {+ :j 1}}
   else :a}}}

{def f1
 {lambda {:a :i :y :x :k}
  {if {< :k {+ {pow 2 :i} 1}}
   then {f1 :a :i {- 1 {* 2 :y :x}} {- :a {* :x :x}} {+ :k 1}}
   else {- :a {/ :x :y}} }}}

Lua

function leftShift(n,p)
    local r = n
    while p>0 do
        r = r * 2
        p = p - 1
    end
    return r
end

-- main

local MAX_IT = 13
local MAX_IT_J = 10
local a1 = 1.0
local a2 = 0.0
local d1 = 3.2

print(" i       d")
for i=2,MAX_IT do
    local a = a1 + (a1 - a2) / d1
    for j=1,MAX_IT_J do
        local x = 0.0
        local y = 0.0
        for k=1,leftShift(1,i) do
            y = 1.0 - 2.0 * y * x
            x = a - x * x
        end
        a = a - x / y
    end
    d = (a1 - a2) / (a - a1)
    print(string.format("%2d    %.8f", i, d))
    d1 = d
    a2 = a1
    a1 = a
end
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

M2000 Interpreter

Using decimal type (26 decimal places) is better for this calculation (this has the same output as FORTRAN). Variable maxitj can be change to lower values when the i get higher value. Although we can't go lower than 2. So here we can start with 13, and lower to 2 for 16th iteration of i

module Feigenbaum_constant_calculation (maxit as integer, c as single){
	locale 1033  // show dot for decimal separator symbol
	single maxitj=13
	integer i, j
	long k
	decimal a1=1, a2, d , d1=3.2, y, x, a
	print "Feigenbaum constant calculation:"
	print
	print format$("{0:-7} {1:-12} {2}","i", "δ","max j")
	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}:a-=x/y}
		d=(a1-a2)/(a-a1)
		print format$("{0::-7} {1:10:-12} {2::-5}",i, d, j-1)
		maxitj-=c
		d1=d:a2=a1:a1= a
	next
}
profiler
Feigenbaum_constant_calculation 18, .7
print timecount
Output:
       i            δ max j
       2 3.2185114220    13
       3 4.3856775986    12 
       4 4.6009492765    12
       5 4.6551304954    11 
       6 4.6661119478    10
       7 4.6685485814    10
       8 4.6690606606     9
       9 4.6691715554     8
      10 4.6691951560     7
      11 4.6692002291     7
      12 4.6692013133     6
      13 4.6692015458     5
      14 4.6692015955     5
      15 4.6692016062     4
      16 4.6692016085     3
      17 4.6692016090     3
      18 4.6692016091     2

Mathematica/Wolfram Language

Translation of: D
maxit = 13;
maxitj = 10;
a1 = 1.0;
a2 = 0.0;
d1 = 3.2;
a = 0.0;
Table[
  a = a1 + (a1 - a2)/d1;
  Do[
   x = 0.0;
   y = 0.0;
   Do[
    y = 1.0 - 2.0 y x;
    x = a - x x;
    ,
    {k, 1, 2^i}
    ];
   a = a - x/y
   ,
   {j, maxitj}
   ];
  d = (a1 - a2)/(a - a1);
  d1 = d;
  a2 = a1;
  a1 = a;
  {i, d}
  ,
  {i, 2, maxit}
] // Grid
Output:
2	3.21851
3	4.38568
4	4.60095
5	4.65513
6	4.66611
7	4.66855
8	4.66906
9	4.66917
10	4.6692
11	4.6692
12	4.6692
13	4.66921

Modula-2

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.

Nim

Translation of: Kotlin
import strformat

iterator feigenbaum(): tuple[n: int; δ: float] =
  ## Yield successive approximations of Feigenbaum constant.

  const
    MaxI = 13
    MaxJ = 10
  var
    a1 = 1.0
    a2 = 0.0
    δ = 3.2

  for i in 2..MaxI:
    var a = a1 + (a1 - a2) / δ
    for j in 1..MaxJ:
      var x, y = 0.0
      for _ in 1..(1 shl i):
        y = 1 - 2 * y * x
        x = a - x * x
      a -= x / y

    δ = (a1 - a2) / (a - a1)
    a2 = a1
    a1 = a
    yield (i, δ)

echo " i         δ"
for n, δ in feigenbaum():
  echo fmt"{n:2d}    {δ:.8f}"
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.66917155
10    4.66919515
11    4.66920026
12    4.66920098
13    4.66920537

Perl

use strict;
use warnings;
use Math::AnyNum 'sqr';

my $a1 = 1.0;
my $a2 = 0.0;
my $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 - sqr($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.60094927653808
 5  4.65513049539198
 6  4.66611194782857
 7  4.66854858144684
 8  4.66906066064827
 9  4.66917155537951
10  4.66919515603002
11  4.66920022908686
12  4.66920131329420
13  4.66920154578091

Phix

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

Python

Translation of: D
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
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

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

Raku

(formerly Perl 6)

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

REXX

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*/
                                                 /*stick a fork in it,  we're all done. */
say left('', 9 + 1 + 11 + 1 + t )"↑"             /*show position of greatest accuracy.  */
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

# 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

RPL

Translation of: Python
Works with: Halcyon Calc version 4.2.7
RPL code Python code
 ≪ { }
  13 10  → maxit maxitj 
  ≪ 3.2 1 0 
     2 maxit 1 + FOR ii 
        DUP2 - 4 PICK / 3 PICK +
        1 maxitj 1 + START
           0 0
           1 2 ii ^ START
              OVER * 2 * 1 SWAP -
              3 PICK ROT SQ - SWAP NEXT
           / - NEXT
        3 PICK ROT - OVER 4 PICK - / 
        5 ROLL OVER + 5 ROLLD
        4 ROLL DROP SWAP ROT NEXT 
     3 DROPN
≫ ≫ 
´FBAUM’ STO 
FBAUM ( -- { δ1..δ13 } )
max_it, max_it_j = 13, 10
d1, a1, a2 = 3.2, 1, 0
for i in range(2, max_it + 1):
   a = a1 + (a1 - a2) / d1
   for j in range(1, max_it_j + 1):
       x = y = 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(d)
   d1, a2, a1 = d, a1, a
clean stack
.
.
Output:
1: { 3.21851142204 4.38567759857 4.60094927654 4.65513049539 4.66611194782 4.66854858152 4.66906066029 4.66917155686 4.6691951528 4.66920033694 4.66920090912 4.66920429563 4.66917851362 }

The above program (limited at 10 iterations) takes 33 minutes and 50 seconds to be executed on a HP-28S.

Ruby

Translation of: C#
def main
    maxIt = 13
    maxItJ = 10
    a1 = 1.0
    a2 = 0.0
    d1 = 3.2
    puts " i       d"
    for i in 2 .. maxIt
        a = a1 + (a1 - a2) / d1
        for j in 1 .. maxItJ
            x = 0.0
            y = 0.0
            for k in 1 .. 1 << i
                y = 1.0 - 2.0 * y * x
                x = a - x * x
            end
            a = a - x / y
        end
        d = (a1 - a2) / (a - a1)
        print "%2d    %.8f\n" % [i, d]
        d1 = d
        a2 = a1
        a1 = a
    end
end

main()
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

Scala

Imperative, ugly

object Feigenbaum1 extends App {
  val (max_it, max_it_j) = (13, 10)
  var (a1, a2, d1, a) = (1.0, 0.0, 3.2, 0.0)

  println(" i       d")
  var i: Int = 2
  while (i <= max_it) {
    a = a1 + (a1 - a2) / d1
    for (_ <- 0 until max_it_j) {
      var (x, y) = (0.0, 0.0)
      for (_ <- 0 until 1 << i) {
        y = 1.0 - 2.0 * y * x
        x = a - x * x
      }
      a -= x / y
    }
    val d: Double = (a1 - a2) / (a - a1)
    printf("%2d    %.8f\n", i, d)
    d1 = d
    a2 = a1
    a1 = a
    i += 1
  }

}

Functional Style, Tail recursive

Output:
Best seen running in your browser either by ScalaFiddle (ES aka JavaScript, non JVM) or Scastie (remote JVM).
object Feigenbaum2 extends App {
  private val (max_it, max_it_j) = (13, 10)

  private def result = {

    @scala.annotation.tailrec
    def outer(i: Int, d1: Double, a2: Double, a1: Double, acc: Seq[Double]): Seq[Double] = {
      @scala.annotation.tailrec
      def center(j: Int, a: Double): Double = {
        @scala.annotation.tailrec
        def inner(k: Int, end: Int, x: Double, y: Double): (Double, Double) =
          if (k < end) inner(k + 1, end, a - x * x, 1.0 - 2.0 * y * x) else (x, y)

        val (x, y) = inner(0, 1 << i, 0.0, 0.0)
        if (j < max_it_j) {
          center(j + 1, a - (x / y))
        } else a
      }

      if (i <= max_it) {
        val a = center(0, a1 + (a1 - a2) / d1)
        val d: Double = (a1 - a2) / (a - a1)

        outer(i + 1, d, a1, a, acc :+ d)
      } else acc
    }

    outer(2, 3.2, 0, 1.0, Seq[Double]()).zipWithIndex
  }

  println(" i     ≈ δ")
  result.foreach { case (δ, i) => println(f"${i + 2}%2d  $δ%.8f") }

}

Sidef

Translation of: Raku
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

Swift

Translation of: C
import Foundation

func feigenbaum(iterations: Int = 13) {
  var a = 0.0
  var a1 = 1.0
  var a2 = 0.0
  var d = 0.0
  var d1 = 3.2

  print(" i       d")

  for i in 2...iterations {
    a = a1 + (a1 - a2) / d1

    for _ in 1...10 {
      var x = 0.0
      var y = 0.0

      for _ in 1...1<<i {
        y = 1.0 - 2.0 * y * x
        x = a - x * x
      }

      a -= x / y
    }

    d = (a1 - a2) / (a - a1)
    d1 = d
    (a1, a2) = (a, a1)

    print(String(format: "%2d    %.8f", i, d))
  }
}

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

Visual Basic .NET

Translation of: C#
Module Module1

    Sub Main()
        Dim maxIt = 13
        Dim maxItJ = 10
        Dim a1 = 1.0
        Dim a2 = 0.0
        Dim d1 = 3.2
        Console.WriteLine(" i       d")
        For i = 2 To maxIt
            Dim a = a1 + (a1 - a2) / d1
            For j = 1 To maxItJ
                Dim x = 0.0
                Dim y = 0.0
                For k = 1 To 1 << i
                    y = 1.0 - 2.0 * y * x
                    x = a - x * x
                Next
                a -= x / y
            Next
            Dim d = (a1 - a2) / (a - a1)
            Console.WriteLine("{0,2:d}    {1:f8}", i, d)
            d1 = d
            a2 = a1
            a1 = a
        Next
    End Sub

End Module
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

V (Vlang)

Translation of: Go
fn feigenbaum() {
    max_it, max_itj := 13, 10
    mut a1, mut a2, mut d1 := 1.0, 0.0, 3.2
    println(" i       d")
    for i := 2; i <= max_it; i++ {
        mut a := a1 + (a1-a2)/d1
        for j := 1; j <= max_itj; j++ {
            mut x, mut y := 0.0, 0.0
            for k := 1; k <= 1<<u32(i); k++ {
                y = 1.0 - 2.0*y*x
                x = a - x*x
            }
            a -= x / y
        }
        d := (a1 - a2) / (a - a1)
        println("${i:2}    ${d:.8f}")
        d1, a2, a1 = d, a1, a
    }
}
 
fn 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

Wren

Translation of: Ring
Library: Wren-fmt
import "./fmt" for Fmt

var feigenbaum = Fn.new {
    var maxIt = 13
    var maxItJ = 10
    var a1 = 1
    var a2 = 0
    var d1 = 3.2
    System.print(" i       d")
    for (i in 2..maxIt) {
        var a = a1 + (a1 - a2)/d1
        for (j in 1..maxItJ) {
            var x = 0
            var y = 0
            for (k in 1..(1<<i)) {
                y = 1 - 2*y*x
                x = a - x*x
            }
            a = a - x/y
        }
        var d = (a1 - a2)/(a - a1)
        System.print("%(Fmt.d(2, i))    %(Fmt.f(0, d, 8))")
        d1 = d
        a2 = a1
        a1 = a
    }
}

feigenbaum.call()
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

XPL0

Translation of: Wren
def  MaxIt = 13, MaxItJ = 10;
real A, A1, A2, D, D1, X, Y;
int  I, J, K;
[A1:= 1.;  A2:= 0.;  D1:= 3.2;
Text(0, " i       d^m^j");
for I:= 2 to MaxIt do
    [A:= A1 + (A1-A2)/D1;
    for J:= 1 to MaxItJ do
        [X:= 0.;  Y:= 0.;
        for K:= 1 to 1<<I do
            [Y:= 1. - 2.*Y*X;
             X:= A - X*X;
            ];
        A:= A - X/Y;
        ];
    D:= (A1-A2) / (A-A1);
    Format(2, 0);  RlOut(0, float(I));
    Format(5, 8);  RlOut(0, D);
    CrLf(0);
    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

zkl

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