Feigenbaum constant calculation: Difference between revisions
Line 1,842: | Line 1,842: | ||
</pre> |
</pre> |
||
=== Version 2 === |
=== Version 2 === |
||
It was observed that variable a already becomes stable in less than 10 iterations. In below program, version 1 is repeated in a more readable form (Original) and maxj was discarded in favor of the dynamic test on a changed variable a. |
It was observed that variable a already becomes stable in less than 10 iterations. In below program, version 1 is repeated in a more readable form (Original) and maxj was discarded in favor of the dynamic test on a changed variable a (Optimized). |
||
<syntaxhighlight lang="rexx"> |
<syntaxhighlight lang="rexx"> |
||
arg n; if n = '' then n = 30; numeric digits n |
arg n; if n = '' then n = 30; numeric digits n |
Revision as of 16:28, 7 July 2024
You are encouraged to solve this task according to the task description, using any language you may know.
- Task
Calculate the Feigenbaum constant.
- See
-
- Details in the Wikipedia article: Feigenbaum constant.
11l
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
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
# 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
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
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
#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#
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++
#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
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
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#
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
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
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
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
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
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
// 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
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
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
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
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
#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)
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
Version 1
/*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
Version 2
It was observed that variable a already becomes stable in less than 10 iterations. In below program, version 1 is repeated in a more readable form (Original) and maxj was discarded in favor of the dynamic test on a changed variable a (Optimized).
arg n; if n = '' then n = 30; numeric digits n
parse version version; say version; glob. = ''
say 'First Fiegenbaum constant, correct to about 11 decimals'
say 'Using algorithm cf RosettaCode'
say
call time('r'); a = Original(); e = format(time('e'),,3)
say 'Original ' a '('e 'seconds)'
call time('r'); a = Optimized(); e = format(time('e'),,3)
say 'Optimized ' a '('e 'seconds)'
call time('r'); a = TrueValue(); e = format(time('e'),,3)
say 'True value' a '('e 'seconds)'
exit
Original:
procedure expose glob.
/* Outer 2 loops with a fixed value */
numeric digits digits()+2
im = 20; jm = 10
a1 = 1; a2 = 0; d1 = 3.2
do i = 2 to im
a = a1 + (a1-a2)/d1
do j = 1 to jm
x = 0; y = 0
do k = 1 to 2**i
y = 1 - 2*x*y; x = a - x*x
end
a = a - x/y
end
d = (a1-a2) / (a-a1)
parse value d a1 a with d1 a2 a1
end
numeric digits digits()-2
return d+0
Optimized:
procedure expose glob.
/* Center loop stops on achieving desired accuracy */
numeric digits digits()+4; numeric fuzz 4
/* Only outer loop maximum */
im = 20
a1 = 1; a2 = 0; d1 = 3.2
do i = 2 to im
a = a1 + (a1-a2)/d1; v = 0
do forever
x = 0; y = 0
do 2**i
y = 1 - 2*x*y; x = a - x*x
end
a = a - x/y
/* Stop second loop when a does not change anymore */
if a = v then
leave
v = a
end
d = (a1-a2) / (a-a1)
parse value d a1 a with d1 a2 a1
end
numeric digits digits()-4
return d+0
TrueValue:
procedure expose glob.
return 4.66920160910299067185320382046620161725818557747576863274565134300413433021131473713868974402394801381716+0
- output:
Still using 30 digits and 10 for the outer loop, this yields:
REXX-ooRexx_5.0.0(MT)_64-bit 6.05 23 Dec 2022 First Fiegenbaum constant, correct to about 11 decimals Using algorithm cf RosettaCode Original 4.66920160910168221375050076911 (97.186 seconds) Optimized 4.66920160910168170137172844325 (21.927 seconds) True value 4.66920160910299067185320382047 (0.000 seconds) Intermediate results: 2 3.218511422038 3 4.385677598568 4 4.600949276538 5 4.655130495392 6 4.666111947829 7 4.668548581447 8 4.669060660648 9 4.669171555380 10 4.669195156030 11 4.669200229087 12 4.669201313294 13 4.669201545781 14 4.669201595537 15 4.669201606198 16 4.669201608481 17 4.669201608970 18 4.669201609074 19 4.669201609097 20 4.669201609102
Apparently you needed 20 iterations to gain 10-11 correct decimals.
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
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
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
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, δ)
}
- 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
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
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)
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
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
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
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
- Programming Tasks
- Solutions by Programming Task
- 11l
- Ada
- ALGOL 68
- AWK
- BASIC
- BASIC256
- Chipmunk Basic
- Just BASIC
- MSX Basic
- True BASIC
- Yabasic
- C
- C sharp
- C++
- D
- Delphi
- SysUtils,StdCtrls
- EasyLang
- F Sharp
- Factor
- Fortran
- FreeBASIC
- Fōrmulæ
- FutureBasic
- Go
- Groovy
- Haskell
- J
- Java
- Jq
- Julia
- Kotlin
- Lambdatalk
- Lua
- M2000 Interpreter
- Mathematica
- Wolfram Language
- Modula-2
- Nim
- Perl
- Phix
- Python
- Racket
- Raku
- REXX
- Ring
- RPL
- Ruby
- Scala
- Sidef
- Swift
- Visual Basic .NET
- V (Vlang)
- Wren
- Wren-fmt
- XPL0
- Zkl