Interesting question regarding Fermat nested radicals

[Rock910]

I'm trying to investigate the behavior of Fermat numbers as nested radicals.
In particular, I'm trying to find the limit as n approaches infinity of the nested radical: sqrt(F_1 + sqrt(F_2 + sqrt(F_3 + ... sqrt(F_n
Does anyone have any idea what I can do?

[madaco]

Are you looking to prove that it converges, efficiently calculate approximations to it, or find a closed form expression for it?

I think I've convinced myself that it converges (but I might have made a mistake in my reasoning).

[PM 2Ring]

My guess is that it converges, or if it diverges it does so extremely slowly.

A brute-force calculation (using mpmath) of 500 terms with 1000 decimal digits of precision yields

3.398589439226157032953566735957045

[Flumble]

It does converge, since F_k = sqrt(F_(k+1)) < sqrt(F_(k+1)+x) so sqrt(F_1+sqrt(F_2+sqrt(... < sqrt(2*sqrt(2*sqrt(... = sqrt(2)*sqrt(sqrt(2*sqrt(... = Π 2^(1/2)^n = 2^(Σ (1/2)^n ) = 2^2 = 4. (these steps are definitely rigorous )

But err, good luck finding out why it's near 3.398. (would be nice if it turns out to be 2+sqrt(2), though it probably won't)

[gmalivuk]

It does converge, since F_k = sqrt(F_(k+1))...

That right there is false, though.

F_k = 2^(2^k)+1
F_k² = 2^(2^k+1)+2^(2^k+1)+1 = F_k+1 + 2^(2^k+1)

You'd need to show that 2^(2^k+1) is less than or equal to the whole rest of the nested radical chain for the rest of your sketch to work at all.

[Flumble]

Oh, right, there's a +1 in the fermat numbers.

[edit]
What if we loosen it to F_k² = F_k+1 + 2^(2k+1) < 2*F_k+1? That should do it, right?
sqrt(F₁+sqrt(F₂+sqrt(... < sqrt(3*sqrt(3*sqrt(... = sqrt(3)*sqrt(sqrt(3*sqrt(... = Π 3^(1/2)^n = 3^(Σ (1/2)^n ) = 3^2 = 9

[>-)]

here's the corresponding OEIS page and it provides a an alternate equivalent formula which looks much easier to work with and also a efficient algorithm

https://oeis.org/A273580

and a lookup in https://isc.carma.newcastle.edu.au/advancedCalc gives no results ;(

[gmalivuk]

Right, dividing everything by sqrt(2) does make it easier to work with as well as quite a bit more obviously convergent, at least to my eyes.

[PM 2Ring]

here's the corresponding OEIS page and it provides a an alternate equivalent formula which looks much easier to work with and also a efficient algorithm

https://oeis.org/A273580

and a lookup in https://isc.carma.newcastle.edu.au/advancedCalc gives no results ;(

Nice! I like that transformation. In hindsight, it's pretty obvious.

Here's some Python 2 / 3 code that implements that algorithm. First, a simple version that just uses core Python. It prints the value of the transformed expression, the value of the nested Fermat number radical down to F₀=3, as well as Rock910's version that goes down to F₁=5.

Code: Select all


from __future__ import print_function, division

# Index of highest Fermat number used
num = 5

# Table of 2**-(2**n)
t = 0.5
tt = [t]
for _ in range(num):
    t *= t
    tt.append(t)

# Initialize tail to phi
x = (1 + 5 ** .5) / 2
for u in range(num, -1, -1):
    print('%d: %.15f' % (u, x))
    x = (x + tt[u] + 1) ** .5

print('%.15f' % x)
print('%.15f' % (x * 2 ** .5))
print('%.15f' % (2 * x * x - 3))


Output:

Code: Select all

5: 1.618033988749895
4: 1.618033988821844
3: 1.618038703990392
2: 1.619242092458812
1: 1.637602544104891
0: 1.699294719613079
1.788657239275619
2.529543326220399
3.398589439226156


This version uses Python's standard decimal module to perform the calculation to high precision. The last couple of digits are generally not correct.

Code: Select all

''' Evaluate an expression composed of nested square roots of Fermat Numbers

    The nth Fermat number is
    F_n = 2 ** (2 ** n) + 1

    We want to evaluate
    sqrt(F_0 + sqrt(F_1 + sqrt(F_2 + sqrt(F_3 + ...

    Written by PM 2Ring 2017.05.24
'''

from __future__ import print_function
from decimal import Decimal as Dec, getcontext

# Index of highest Fermat number used
num = 7

ctx = getcontext()

# Precision in bits
bin_prec = 1 << (num + 1)
# 10 bits ~= 3 decimal digits
prec = bin_prec * 3 // 10 
ctx.prec = prec
print('precision:', prec)

# Table of 2 ** -(2 ** n)
t = Dec('0.5')
tt = [t]
for _ in range(num):
    t *= t
    tt.append(t)

# Initialize tail to phi
x = (1 + Dec(5).sqrt()) / 2
phi = x
for u in range(num, -1, -1):
    #print(u, x)
    x = (1 + tt[u] + x).sqrt()

# Including F(0)
print(x * Dec(2).sqrt())
# Excluding F(0)
print(2 * x * x - 3)
 

Output

Code: Select all

precision: 76
2.529543326220398430310379128859753335193537124459383417865718711396730946540
3.398589439226157032953566735957045123597355357178507387913820594572190912451


Here's 1000 decimal places of Rock910's expression, calculated using the above code with `num = 11`, which gives a precision of 1228.

Code: Select all

3.
39858 94392 26157 03295 35667 35957 04512 35973 55357 17850
73879 13820 59457 21909 12452 86035 50704 19108 74685 45221
20553 84825 00341 70795 42003 60898 73353 66485 52656 77555
17031 27277 75558 13377 49466 84627 98582 01898 63540 28813
04640 56404 70497 09601 72534 52481 10905 73507 75176 68859

52141 17296 89137 26628 06758 98570 21096 10495 26923 65557
18787 57000 15248 70654 25121 19527 72270 12272 91335 35442
08150 39744 36410 83718 11256 07796 12708 79908 15955 29118
47332 07196 74691 91753 62310 07334 61352 15807 99594 64470
04004 54675 59610 49453 42221 96066 61767 64036 83593 45264

36421 06461 39582 83661 71300 11890 48054 59837 63233 59485
25517 08786 62460 78720 38138 60574 24386 92971 43015 13162
98979 87759 15908 84052 94382 52459 85625 26584 15783 28355
58871 46539 89882 42987 55074 45634 64185 12101 49133 28125
05845 36146 29517 83397 33818 14499 99070 94856 32213 75004

64252 52429 14722 73682 51833 43784 64977 11503 27185 81786
12754 43385 40027 61941 34466 84334 86702 78684 58962 11344
83980 63000 46481 03370 04871 00268 56609 92688 05608 70103
97676 58441 83671 29341 42943 47947 51776 56797 32848 95745
01635 39228 50501 27379 26864 51554 43483 14136 83533 93702


As explained in the OEIS link, this algorithm uses the golden ratio phi as an estimate of the tail of the nested radical. There's not much point actually using a better estimate because the algorithm already converges very quickly. However, it's not that hard to find better estimates.

Here's some code that uses the 3rd-party mpmath module to implement the algorithm. It shows that `phi + 2 ** -(2 ** (n+1)) / (2 * phi)` is a better estimate.

Code: Select all


from mpmath import mp

num = 11

mp.prec = 1 << (num + 1)
print('dps:', mp.dps)

# Table of 2 ** -(2 ** n)
t = mp.mpf('0.5')
tt = [t]
for _ in range(num):
    t *= t
    tt.append(t)

# Number of digits to print
pdps = 20

phi = +mp.phi
x = phi
print('u, x-phi, d=x-z, F(u+2)/d')
for u in range(num, -1, -1):
    if u < num - 1:
        z = phi + tt[u+1] / (2 * phi)
        d = x - z
        c = tt[u+2] / d
        print(u, mp.nstr(x - phi, pdps), mp.nstr(d, pdps), mp.nstr(c, pdps))
    x = mp.sqrt(x + tt[u] + 1)


Output:

Code: Select all

dps: 1232
u, x-phi, d=x-z, tt[u+2]/d
9 1.7189640900454060261e-309 2.041742526503038669e-618 15.155417527999327029
8 2.3047540357797355452e-155 3.6704265230509527055e-310 15.155417527999327029
7 2.6687228498107178181e-78 4.921237384203419835e-156 15.155417527999327029
6 9.0811932798958323201e-40 5.6984035834969891644e-79 15.155417527999327029
5 1.6751844831812945028e-20 1.9390662590631130671e-40 15.155417527999327029
4 7.1948625703880165185e-11 3.5769459021312081437e-21 15.155417528673238374
3 4.7152404967159728467e-6 1.5362820901067768655e-11 15.155461692434825036
2 0.0012081037089170574727 1.0060746399190973412e-6 15.166657081950721751
1 0.019568555354996283934 0.00025499320656206992759 15.319035564381392928
0 0.081260730863183668272 0.0040064822694468122466 15.59971960355875326


Clearly, we could improve the estimate even further using `tt[u+2]/15.155417527999327029`.

[PM 2Ring]

I've done some more work on estimating the tail term of the transformed nested radical. To estimate the tail of the i'th term we can use

Here are the a_j

Spoiler:

0: phi,
1: -2 + 2 * phi,
2: 14 - 8 * phi,
3: -72 + 44 * phi,
4: 410 - 250 * phi,
5: -4340 + 2680 * phi,
6: 44712 - 27644 * phi,
7: -432696 + 267432 * phi,
8: 4434798 - 2740596 * phi,
9: -47623904 + 29432988 * phi,

As decimals:
0: 1.61803398874989
1: 1.23606797749979
2: 1.05572809000084
3: -0.806504495004627
4: 5.49150281252629
5: -3.66891015028181
6: -16.9315850020932
7: 18.065679361879
8: 420.522567993179
9: -329.025532209932

This tail term is much better than phi for large i, and for small i we can use it to approximate the value of the expression without nested radicals.

Setting i = 0 we get

Spoiler:

x = 2 * (phi +
(-2 + 2 * phi) / 16 +
(14 - 8 * phi) / 16 ** 2 +
(-72 + 44 * phi) / 16 ** 3 +
(410 - 250 * phi) / 16 ** 4 +
(-4340 + 2680 * phi) / 16 ** 5 +
(44712 - 27644 * phi) / 16 ** 6 +
(-432696 + 267432 * phi) / 16 ** 7 +
(4434798 - 2740596 * phi) / 16 ** 8 +
(-47623904 + 29432988 * phi) / 16 ** 9)

= 3.398589441285548, which is accurate to 9 significant figures.

Setting i = 1 we get

Spoiler:

(5 + 4 * phi +
(-2 + 2 * phi) / 4 ** 2 +
(14 - 8 * phi) / 4 ** 5 +
(-72 + 44 * phi) / 4 ** 8 +
(410 - 250 * phi) / 4 ** 11 +
(-4340 + 2680 * phi) / 4 ** 14 +
(44712 - 27644 * phi) / 4 ** 17 +
(-432696 + 267432 * phi) / 4 ** 20 +
(4434798 - 2740596 * phi) / 4 ** 23 +
(-47623904 + 29432988 * phi) / 4 ** 26) ** 0.5

= 3.3985894392261576, which is accurate to 15 significant figures, approximately the limit of the IEEE 754 double-precision floating-point format.

Here's some mpmath Python code. The first section of output shows the differences between the true tail term and the successive estimates.Then it shows the result of evaluating several terms of the nested radical, first using phi for the tail, and then using the estimate given above.

Code: Select all


from mpmath import mp
sqrt = mp.sqrt

num = 13

mp.prec = 1 << (num + 2)
print('dps:', mp.dps)

def show(*args, pdps=50):
    print(*[mp.nstr(u, pdps) for u in args])

# Some constants
two = mp.mpf(2)
four = mp.mpf(4)
phi = +mp.phi

def tt(n):
    return two ** -(two ** n)

# Coefficients used to estimate the tail of the nested radical
aa = [
    phi,
    -2 + 2 * phi,
    14 - 8 * phi,
    -72 + 44 * phi,
    410 - 250 * phi,
    -4340 + 2680 * phi,
    44712 - 27644 * phi,
    -432696 + 267432 * phi,
    4434798 - 2740596 * phi,
    -47623904 + 29432988 * phi,
]

#for i, a in enumerate(aa):
    #show(i, a, pdps=15)
#print()

x = phi
for i in range(num, -1, -1):
    z = 0
    dd = []
    k = two**i + 1
    for j, a in enumerate(aa):
        z += a * four ** -(j * k)
        dd.append(x - z)
    print(i, end=': ')
    show(*dd, pdps=5)
    x = sqrt(x + tt(i) + 1)
print('\n*', end=' ')
show(2 * x * x - 3, pdps=120)

print('\nplain')
for i in range(7):
    x = phi
    for j in range(i, 0, -1):
        x = sqrt(x + tt(j) + 1)
    show(i, 2 * x, pdps=2+2**i)

print('\nenhanced')
for i in range(5):
    x = 0
    k = two**i + 1
    for j, a in enumerate(aa):
        x += a * four ** -(j * k)
    for j in range(i, 0, -1):
        x = sqrt(x + tt(j) + 1)
    show(i, 2 * x, pdps=2**(3+i))


Output

Code: Select all

dps: 9863
13: 0.0 -2.5974e-4933 -2.5974e-4933 -2.5974e-4933 -2.5974e-4933 -2.5974e-4933 -2.5974e-4933 -2.5974e-4933 -2.5974e-4933 -2.5974e-4933
12: 2.8331e-2467 -2.4803e-4934 -8.0263e-4934 -8.0263e-4934 -8.0263e-4934 -8.0263e-4934 -8.0263e-4934 -8.0263e-4934 -8.0263e-4934 -8.0263e-4934
11: 2.9588e-1234 6.0493e-2468 -1.1062e-3701 -6.7724e-4935 -2.4803e-4934 -2.4803e-4934 -2.4803e-4934 -2.4803e-4934 -2.4803e-4934 -2.4803e-4934
10: 9.5621e-618 6.3179e-1235 -3.7336e-1852 1.9666e-2468 -1.0164e-3085 -3.6287e-3702 2.9951e-4319 -2.2711e-4935 -7.6644e-4935 -7.6644e-4935
9: 1.719e-309 2.0417e-618 -2.1691e-927 2.0539e-1235 -1.9084e-1544 -1.2247e-1852 1.8173e-2161 5.8828e-2469 -6.401e-2778 -3.0347e-3086
8: 2.3048e-155 3.6704e-310 -5.2282e-465 6.6377e-619 -8.2689e-774 -7.1153e-928 1.4156e-1082 6.1439e-1236 -8.9633e-1391 -5.6976e-1545
7: 2.6687e-78 4.9212e-156 -8.1169e-234 1.1933e-310 -1.7212e-388 -1.715e-465 3.9508e-543 1.9855e-619 -3.3541e-697 -2.4688e-774
6: 9.0812e-40 5.6984e-79 -3.1982e-118 1.5999e-156 -7.8531e-196 -2.6626e-234 2.0872e-273 3.5694e-311 -2.0518e-350 -5.139e-389
5: 1.6752e-20 1.9391e-40 -2.0076e-60 1.8526e-79 -1.6774e-99 -1.0491e-118 1.517e-138 4.7858e-157 -5.0747e-177 -2.3446e-196
4: 7.1949e-11 3.5769e-21 -1.5905e-31 6.3039e-41 -2.4515e-51 -6.5853e-61 4.0899e-71 5.5415e-80 -2.5238e-90 -5.0081e-100
3: 4.7152e-6 1.5363e-11 -4.4769e-17 1.1629e-21 -2.9638e-27 -5.2174e-32 2.1238e-37 1.8857e-41 -5.6283e-47 -7.3193e-52
2: 0.0012081 1.0061e-6 -7.4612e-10 4.9912e-12 -3.2733e-15 -1.467e-17 1.565e-20 3.4758e-22 -2.6667e-25 -8.8396e-28
1: 0.019569 0.00025499 -2.7529e-6 3.2366e-7 -3.6577e-9 -2.408e-10 5.5824e-12 1.4748e-12 -1.9226e-14 -9.6175e-16
0: 0.081261 0.0040065 -0.00011746 7.9445e-5 -4.3488e-6 -8.4981e-7 1.5939e-7 9.2093e-8 -5.8176e-9 -1.0297e-9

* 3.3985894392261570329535667359570451235973553571785073879138205945721909124528603550704191087468545221205538482500341708

plain
0 3.24
1 3.387
2 3.39837
3 3.398589178
4 3.39858943922492349
5 3.398589439226157032953477984168113
6 3.39858943922615703295356673595704512359735387042231065900885926106

enhanced
0 3.3985894
1 3.398589439226158
2 3.3985894392261570329535667361159
3 3.398589439226157032953566735957045123597355357178507428522643301
4 3.3985894392261570329535667359570451235973553571785073879138205945721909124528603550704191087468545221291401062910876379386856078