It was in January of 2017 that I made my first and only post about one of the two Feigenbaum constants 4.66920160910299067185320382... The YouTube video that I included in that post is still active and I've embedded it in this post as well. Watching that will explain the significance of this constant. Here is a link to another video by Veritasium that covers the same ground.
A257809 Lexicographically largest strictly increasing sequence of primes for which the continued square root map produces Feigenbaum's constant delta = 4.6692016... (A006890).
Computing a Prime Sequence using Feigenbaum's Constant in SageMath
Here is the complete SageMath code for generating a sequence of primes dictated by deeply nested square roots, incorporating Feigenbaum's constant. A highly precise value for $\delta$ is provided to ensure the sequence evaluates correctly, as nested square root calculations can be sensitive to floating-point precision.
The SageMath Code
# Force 500 bits of precision (approx 150 decimal digits)
R = RealField(500)
# Feigenbaum's constant delta extended to 100 decimal digits
delta = R("4.6692016091029906718532038204662016172581855774757686327456513430041343302113147371386897440239480138")
def CSR(v, s=0):
# Ensure our running sum uses the 500-bit precision field
s = R(s)
for x in reversed(v):
s = R(x + s).sqrt()
return s
a = [13]
for n in range(50):
print(a[-1], end=", ")
i = prime_pi(a[-1]) + 1
while True:
next_primes = [nth_prime(i + j) for j in range(1, 10)]
v_test = a + next_primes
if CSR(v_test) >= delta:
a.append(nth_prime(i))
break
i += 1
The Mathematical Context
Feigenbaum's constant $\delta \approx 4.6692016$ is a universal mathematical constant that arises in chaos theory and non-linear dynamics. Specifically, it describes the rate at which period-doubling bifurcations occur as a system approaches chaotic behavior.
While the constant was originally discovered in the context of the logistic map and population dynamics, in this code, it acts purely as a strict numeric threshold. The algorithm tests candidate prime numbers by nesting them inside a Continued Square Root (CSR) function. The sequence only accepts a new prime when the deeply nested square root expression of the sequence exceeds this chaotic boundary.
Because the sequence uses 9 lookahead primes and deeply nested roots, SageMath's RR (which defaults to 53 bits of precision) is perfectly suited to handle the threshold check rapidly without getting bogged down by the computational weight of exact symbolic algebra.
Step-by-Step Manual Trace
To understand the progression, we have to look at how Feigenbaum's constant dictates the threshold at each layer. Squaring $\delta$ gives us our initial internal target: $\delta^2 \approx 21.80144$.
Iteration 1
- Current sequence:
[13] - The Condition: To append a new prime, the algorithm requires the nested root of the sequence plus a 9-prime lookahead to cross the threshold:
$$\sqrt{13 + \sqrt{p_1 + \dots}} \ge 4.6692016$$ - Unpacking the Math: Squaring both sides and subtracting 13 strips away the outer layer:
$$\sqrt{p_1 + \dots} \ge 8.80144$$
Squaring once more gives the approximate threshold the lookahead primes must reach:
$$p_1 + \dots \ge 77.465$$ - The Search:
- If the script tests the 18th prime (61), its lookahead sequence starts with 67. The continued square root of this lookahead evaluates to $\approx 8.72$, which falls short of the 8.80144 requirement.
- When the script tests the 19th prime (67), the lookahead sequence starts with 71. The continued square root of this sequence evaluates to $\approx 8.95$, successfully crossing the threshold.
- Result: The 19th prime is appended.
- New Sequence:
[13, 67]
Iteration 2
- Current sequence:
[13, 67] - The Condition: The nested expression is now one layer deeper:
$$\sqrt{13 + \sqrt{67 + \sqrt{p_1 + \dots}}} \ge 4.6692016$$ - Unpacking the Math: We already know from Iteration 1 that to satisfy the outer layer, the inner contents must exceed 8.80144:
$$\sqrt{67 + \sqrt{p_1 + \dots}} \ge 8.80144$$
Squaring this and subtracting 67 strips the second layer:
$$\sqrt{p_1 + \dots} \ge 10.4654$$
Squaring again shows what the next lookahead sequence must achieve:
$$p_1 + \dots \ge 109.52$$ - The Search:
- Testing the 24th prime (89) uses a lookahead starting at 97, yielding a root of $\approx 10.37$ (too small).
- Testing the 25th prime (97) uses a lookahead starting at 101, yielding a root of $\approx 10.56$.
- Result: The 25th prime is appended.
- New Sequence:
[13, 67, 97]
Iteration 3
- Current sequence:
[13, 67, 97] - The Condition: The expression adds another layer:
$$\sqrt{13 + \sqrt{67 + \sqrt{97 + \sqrt{p_1 + \dots}}}} \ge 4.6692016$$ - Unpacking the Math: Relying on the unpacked threshold from Iteration 2:
$$\sqrt{97 + \sqrt{p_1 + \dots}} \ge 10.4654$$
Squaring and subtracting 97 gives the new inner requirement:
$$\sqrt{p_1 + \dots} \ge 12.525$$
Squaring this dictates the next lookahead target:
$$p_1 + \dots \ge 156.88$$ - The Search:
- Testing the 33rd prime (137) yields a lookahead root of $\approx 12.3$.
- Testing the 34th prime (139) uses a lookahead starting at 149, yielding a root of $\approx 12.7$.
- Result: The 34th prime is appended.
- New Sequence:
[13, 67, 97, 139]
Asymptotic Growth Rate Script
Because of the nested square roots, the values required to maintain the threshold grow extremely quickly. The following script calculates the sequence, analyzes its growth rate, and visualizes the results on a logarithmic scale.
# Feigenbaum's constant delta
delta = 4.66920160910299067185320382
def CSR(v, s=0):
for x in reversed(v):
s = RR(x + s).sqrt()
return s
a = [13]
growth_rates = []
# Generate the first 15 terms (kept smaller as it grows massively)
for n in range(15):
i = prime_pi(a[-1]) + 1
while True:
next_primes = [nth_prime(i + j) for j in range(1, 10)]
v_test = a + next_primes
if CSR(v_test) >= delta:
new_prime = nth_prime(i)
# Calculate the ratio between the new prime and the square of the previous
if len(a) > 0:
growth_rates.append(RR(new_prime) / RR(a[-1]^2))
a.append(new_prime)
break
i += 1
print("Sequence:", a)
print("Growth ratios (a_n / a_{n-1}^2):", [round(g, 4) for g in growth_rates])
# Plotting the sequence on a logarithmic scale
list_plot(a, scale='semilogy', axes_labels=['Index (n)', 'Prime Value (Log Scale)'],
title="Asymptotic Growth of the Sequence", color='blue', size=30)
Asymptotic Growth of the Sequence
The Asymptotic Analysis
When dealing with deeply nested square roots, the dominant mathematical constraint acting on the sequence is the squaring required to "unpeel" each layer.
If we approximate the relationship between two adjacent terms $a_n$ and $a_{n+1}$ to maintain a constant threshold $k$, the equation simplifies to:
$$\sqrt{a_n + \sqrt{a_{n+1} + \dots}} \approx k$$
Squaring both sides and isolating the inner root gives:
$$\sqrt{a_{n+1} + \dots} \approx k^2 - a_n$$
To keep the right side of the equation positive and stable, the next term in the sequence, $a_{n+1}$, must scale quadratically with respect to the current term $a_n$. This means the relationship is roughly:
$$a_{n+1} \approx (a_n)^2$$
Because each term is the square of the previous term, the sequence exhibits double exponential growth. The general closed-form asymptotic behavior for a sequence of this nature is:
$$a_n \sim C^{2^n}$$
where $C$ is a constant defined by the initial conditions of the system. In standard exponential growth, a variable increases by a fixed percentage. In double exponential growth, the number of digits in the number increases exponentially, which is why generating terms of this specific sequence requires significant computational power and highly optimized floating-point arithmetic.
