Super Golden Ratio

[Thomas Klemm]

Recently I stumbled upon this video by Micheal Penn: What is the super-golden ratio??

It is defined by this equation:

\(
\begin{align}
\psi^3=\psi^2+1
\end{align}
\)

He solves it algebraically to find:

\(
\begin{align}
\psi=\frac{1}{3}\left(1+\sqrt[3]{\frac{1}{2}\left(29-3\sqrt{93}\right)}+\sqrt[3]{\frac{1}{2}\left(29+3\sqrt{93}\right)}\right)
\end{align}
\)

We can calculate a numerical approximation using Bernoulli's Method with the following program for the HP-42S:

Code:

00 { 4-Byte Prgm }
01 ENTER
02 RCL+ ST T
03 END

Initialise the stack with:

CLST
1

And then hit repeatedly the R/S key:

R/S
R/S
R/S
…

It produces the following sequence:

1 1 1 2 3 4 6 9 13 19 28 …

This is also known as A000930: Narayana's cows sequence.

The quotients of consecutive elements converge to the root and thus to \(\psi\):

Code:

           1 ÷            1 =  1.00000000000
           2 ÷            1 =  2.00000000000
           3 ÷            2 =  1.50000000000
           4 ÷            3 =  1.33333333333
           6 ÷            4 =  1.50000000000
           9 ÷            6 =  1.50000000000
          13 ÷            9 =  1.44444444444
          19 ÷           13 =  1.46153846154
          28 ÷           19 =  1.47368421053
          41 ÷           28 =  1.46428571429
          60 ÷           41 =  1.46341463415
          88 ÷           60 =  1.46666666667
         129 ÷           88 =  1.46590909091
         189 ÷          129 =  1.46511627907
         277 ÷          189 =  1.46560846561
         406 ÷          277 =  1.46570397112
         595 ÷          406 =  1.46551724138
         872 ÷          595 =  1.46554621849
        1278 ÷          872 =  1.46559633028
        1873 ÷         1278 =  1.46557120501
        2745 ÷         1873 =  1.46556326749
        4023 ÷         2745 =  1.46557377049
        5896 ÷         4023 =  1.46557295551
        8641 ÷         5896 =  1.46556987788
       12664 ÷         8641 =  1.46557111445
       18560 ÷        12664 =  1.46557169931
       27201 ÷        18560 =  1.46557112069
       39865 ÷        27201 =  1.46557111871
       58425 ÷        39865 =  1.46557130315
       85626 ÷        58425 =  1.46557124519
      125491 ÷        85626 =  1.46557120501
      183916 ÷       125491 =  1.46557123618
      269542 ÷       183916 =  1.46557123904
      395033 ÷       269542 =  1.46557122823
      578949 ÷       395033 =  1.46557123076
      848491 ÷       578949 =  1.46557123339
     1243524 ÷       848491 =  1.46557123175
     1822473 ÷      1243524 =  1.46557123144
     2670964 ÷      1822473 =  1.46557123206
     3914488 ÷      2670964 =  1.46557123196
     5736961 ÷      3914488 =  1.46557123179
     8407925 ÷      5736961 =  1.46557123188
    12322413 ÷      8407925 =  1.46557123190
    18059374 ÷     12322413 =  1.46557123187
    26467299 ÷     18059374 =  1.46557123187
    38789712 ÷     26467299 =  1.46557123188
    56849086 ÷     38789712 =  1.46557123188
    83316385 ÷     56849086 =  1.46557123188
   122106097 ÷     83316385 =  1.46557123188
   178955183 ÷    122106097 =  1.46557123188

But this is the power iteration to calculate the greatest (in absolute value) eigenvalue of a diagonalizable matrix.

To see this enter the following matrix with:

3 ENTER NEWMAT

\(
\begin{bmatrix}
1 & 0 & 1 \\
1 & 0 & 0 \\
0 & 1 & 0 \\
\end{bmatrix}
\)

Fill the stack using:

ENTER
ENTER
ENTER

Enter the initial vector:

\(
\begin{bmatrix}
1 \\
1 \\
1 \\
\end{bmatrix}
\)

3 ENTER 1 NEWMAT
1 +

And now repeatedly hit the multiply key:

*
*
*
…

If you think you had enough, calculate the quotient of two consecutive elements, e.g. for:

\(
\begin{bmatrix}
58425 \\
39865 \\
27201 \\
\end{bmatrix}
\)

58425
39865
÷

1.46557130315 (1.46557123188)

This is a "Good place to stop".

---

References

You can find a description in Computational analysis with the HP-25 pocket calculator (Peter Henrici):

• Bernoulli's Method for Single Dominant Zero (p. 80)
  

[Thomas Klemm]

If you want some more do the same with the Plastic number or the Tribonacci numbers.

• The Plastic Ratio - Numberphile
• Phi and the TRIBONACCI monster
  

Or then use one of the many algorithms to accelerate the convergence of the original sequence.

[Albert Chan]

“Rabbit Math” – A Formula for the Magical Fibonacci Sequence

We use the same trick, for "Super" Fibonacci numbers generation function

x/(1-x-x^3) = x + x^2 + x^3 + 2*x^4 + 3*x^5 + 4*x^6 + 6*x^7 + 9*x^8 + ...

Binet-like direct formula for "Super" Fibonacci numbers, has this form:

S(n) = [k1, k2, k3] * [r1, r2, r3] .^ n

[r1, r2, r3] are reciprocal roots of 1-x-x^3 = 0 (or, roots of x^3 = x^2+1)
[k1, k2, k3] are constants to be solved, so that S(1) = S(2) = S(3) = 1

XCAS> R := 1 ./ proot([-1,0,-1,1]);
→ [-0.232786+0.792552*i, -0.232786-0.792552*i, 1.46557]

XCAS> K := proot([-1,0,3/31,1/31]);
→ [-0.208619-0.183825*i, -0.208619+0.183825*i, 0.417238]

XCAS> [K*R.^29, K*R.^30, K*R.^31]
→ [27201.0, 39865.0, 58425.0]

Both cubics are in depressed form, easy to get (K,R) exact form

w = exp(i*2*pi/3) :;
kb := (1/62 + sqrt(93)*3/1922)^(1/3) :;
ka := +1/(31*kb) :;
rb := (1/2 + sqrt(93)/18)^(1/3) :;
ra := -1/(3*rb) :;

K := [ka*w+kb/w, ka/w+kb*w, ka+kb] :;
R := 1 ./ [ra*w+rb/w, ra/w+rb*w, ra+rb] :;

approx(K,R) matched proot values.
We just need to confirm K coefficients produce S(1) = S(2) = S(3) = 1, exactly.

XCAS> S(n) := K * R.^n
XCAS> simplify([S(1), S(2), S(3)])

After 31 seconds, we get back [1, 1, 1]

[Gerson W. Barbosa]

(08-07-2022 07:55 AM)Thomas Klemm Wrote:  We can calculate a numerical approximation using Bernoulli's Method with the following program for the HP-42S:

Code:

00 { 4-Byte Prgm }
01 ENTER
02 RCL+ ST T
03 END

Initialise the stack with:

CLST
1

And then hit repeatedly the R/S key:

R/S
R/S
R/S
…

This will do for the 12 digits in the display:

67
ENTER
X↑2
1/X
+
11
1/X
Y↑X

Or, for more digits,

(67 + 1/(67^2 + 2/(67 + 3/(2×67^2 + 11/(3×67)))))^(1/11)
= 1.46557123187676802665673122475

The latter should take up more steps than the exact expression, though.

P.S.:

300766 ENTER 33 1/X Y↑X is shorter for 12 digits and the approximation is slightly better.

[Thomas Klemm]

(08-08-2022 04:07 PM)Albert Chan Wrote:  S(n) = [k1, k2, k3] * [r1, r2, r3] .^ n

[r1, r2, r3] are reciprocal roots of 1-x-x^3 = 0 (or, roots of x^3 = x^2+1)
[k1, k2, k3] are constants to be solved, so that S(1) = S(2) = S(3) = 1

We assume that \(\alpha\), \(\beta\) and \(\gamma\) are the roots of \(x^3=x^2+1\).
Thus we have:

\(
\begin{align}
\alpha + \beta + \gamma &= 1 \\
\alpha \cdot \beta \cdot \gamma &= 1 \\
\end{align}
\)

For the direct formula \(S(n) = u \cdot \alpha^n + v \cdot \beta^n + w \cdot \gamma^n\) the constants \(u\), \(v\) and \(w\) are solved so that: \(S(1) = S(2) = S(3) = 1\)

This leads to:

\(
\begin{bmatrix}
\alpha & \beta & \gamma \\
\alpha^2 & \beta^2 & \gamma^2 \\
\alpha^3 & \beta^3 & \gamma^3 \\
\end{bmatrix}
\begin{bmatrix}
u \\
v \\
w \\
\end{bmatrix}
=
\begin{bmatrix}
1 \\
1 \\
1 \\
\end{bmatrix}
\)

This can be solved for \(u\) to get:

\(
\begin{align}
u
&= \frac{(1-\beta)(1-\gamma)}{\alpha(\alpha-\beta)(\alpha-\gamma)} \\
&= \frac{1 - (\beta+\gamma) + \beta \gamma}{a[\alpha^2 - (\beta+\gamma)\alpha + \beta \gamma]} \\
&= \frac{1 - (1-\alpha) + \frac{1}{\alpha}}{\alpha[\alpha^2 - (1-\alpha)\alpha + \frac{1}{\alpha}]} \\
&= \frac{\alpha + \frac{1}{\alpha}}{\alpha[2\alpha^2 - \alpha + \frac{1}{\alpha}]} \\
&= \frac{1+ \frac{1}{\alpha^2}}{2\alpha^2 - \alpha + \frac{1}{\alpha}} \\
&= \frac{1+ \frac{1}{\alpha^2}}{2\alpha^2 - 2\alpha + \alpha +\frac{1}{\alpha}} \\
&= \frac{\alpha}{3\alpha^2 - 2\alpha} \\
&= \frac{1}{3\alpha - 2}
\end{align}
\)

And similarly for \(v\) and \(w\).
We consider \(\alpha = \psi\), i.e. the real solution which is dominant.

The following formula can be used to calculate \(\psi\):

\(
\begin{align}
\psi ={\frac {2}{3}}\cosh {\left({\tfrac {\cosh ^{-1}\left({\frac {29}{2}}\right)}{3}}\right)}+{\frac {1}{3}}
\end{align}
\)

This allows to calculate both \(\alpha\) in R00 and \(\frac{1}{u}\) in R01:

Code:

00 { 30-Byte Prgm }
01 29
02 2
03 ÷
04 ACOSH
05 3
06 ÷
07 COSH
08 2
09 ×
10 1
11 +
12 RCL ST X
13 3
14 ÷
15 STO 00
16 X<>Y
17 2
18 -
19 STO 01
20 END

R00: 1.465571231876768026656731225219939
R01: 2.396713695630304079970193675659816

We cheat a little bit and use only the dominant root and round the result to the next integer:

Code:

00 { 12-Byte Prgm }
01 RCL 00
02 X<>Y
03 Y↑X
04 RCL÷ 01
05 0.5
06 +
07 IP
08 END

Examples

1 R/S
1

2 R/S
1

3 R/S
1

4 R/S
2

5 R/S
3

6 R/S
4

7 R/S
6

29 R/S
27201

30 R/S
39865

31 R/S
58425

---

(08-08-2022 04:07 PM)Albert Chan Wrote:  XCAS> K := proot([-1,0,3/31,1/31]);
→ [-0.208619-0.183825*i, -0.208619+0.183825*i, 0.417238]

We can verify that indeed: \(31u^3=3u+1\)

\(
\begin{align}
\frac{3}{u^2} + \frac{1}{u^3}
&= 3(3\alpha-2)^2 + (3\alpha-2)^3 \\
&= (3\alpha-2)^2(3\alpha+1) \\
&= (9\alpha^2-12\alpha+4)(3\alpha+1) \\
&= 27\alpha^3 - 27\alpha^2 + 4 \\
&= 27(\alpha^3 - \alpha^2) + 4 \\
&= 27 + 4 \\
&= 31
\end{align}
\)

How did you come up with this equation?

[Albert Chan]

(08-12-2022 05:46 PM)Thomas Klemm Wrote:  We assume that \(\alpha\), \(\beta\) and \(\gamma\) are the roots of \(x^3=x^2+1\).
Thus we have:

\(
\begin{align}
\alpha + \beta + \gamma &= 1 \\
\alpha \cdot \beta \cdot \gamma &= 1 \\
\end{align}
\)

For the direct formula \(S(n) = u \cdot \alpha^n + v \cdot \beta^n + w \cdot \gamma^n\) the constants \(u\), \(v\) and \(w\) are solved so that: \(S(1) = S(2) = S(3) = 1\)

This leads to:

\(
\begin{bmatrix}
\alpha & \beta & \gamma \\
\alpha^2 & \beta^2 & \gamma^2 \\
\alpha^3 & \beta^3 & \gamma^3 \\
\end{bmatrix}
\begin{bmatrix}
u \\
v \\
w \\
\end{bmatrix}
=
\begin{bmatrix}
1 \\
1 \\
1 \\
\end{bmatrix}
\)

This can be solved for \(u\) to get ...

Simpler approach is do partial fraction decomposition of generation function.
S(n) formula is coefficient of x^n, with RHS geometric series in normalized form.

\(\displaystyle
\frac{x}{(1 - α x)(1- β x) (1 - γ x)} =
\frac{u}{1 - α x} +
\frac{v}{1- β x} +
\frac{w}{1 - γ x}
\)

Mulitply both side by (1 - α x), then let x = 1/α, we solved u

\(\displaystyle u = \frac{1/α}{(1 - β/α) (1 - γ/α)}
= \frac{α}{(α-β)(α-γ)}\)

\(\displaystyle \frac{1}{u} = α - (β+γ) + \frac{βγ}{α} = α - (1-α) + (α-1) = (3α - 2)\)

\(\displaystyle \frac{1}{u} + \frac{1}{v} + \frac{1}{w}
= (3α\!-\!2) + (3β\!-\!2) + (3γ\!-\!2) = 3 - 6 = -3 \)

If we let x = 0, we get 0 = u + v + w \(\displaystyle \;⇒ \frac{1}{uv} + \frac{1}{uw} + \frac{1}{vw} = 0\)

\(\displaystyle \frac{1}{uvw}
= (3α\!-\!2)(3β\!-\!2)(3γ\!-\!2) = 27(αβγ) - 18(αβ\!+\!αγ\!+\!βγ) + 12(α\!+\!β\!+\!γ) - 8
= 27 + 12 - 8 = 31
\)

\(\displaystyle ⇒ (x-\frac{1}{u})(x-\frac{1}{v})(x-\frac{1}{w}) = x^3 + 3x^2 - 31\)

[Albert Chan]

Another way to get cubic with roots = (roots of x^3-x^2-1) * 3 - 2

CAS> subst(x^3 - x^2 - 1, x = (y+2)/3)

((y+2)/3)^3 - ((y+2)/3)^2 - 1

CAS> numer(Ans) // y has roots 1/u, 1/v, 1/w

y^3 + 3*y^2 - 31

[Thomas Klemm]

Thanks a lot for your explanations.

[Albert Chan]

(08-09-2022 10:45 AM)Gerson W. Barbosa Wrote:  300766 ENTER 33 1/X Y↑X is shorter for 12 digits and the approximation is slightly better.

Let (α,β,γ) be roots of x^3 = x^2 + 1
We can build cubic with roots (-α^n, -β^n, -γ^n)

n=1: x^3 + x^2 + 1
n=2: x^3 + x^2 - 2x + 1
n=3: x^3 + 4x^2 + 3x + 1

x^2 term coefficient = (α^+n + β^+n + γ^+n)
x^1 term coefficient = (α^−n + β^−n + γ^−n), since αβγ = 1

x^+n = x^+(n-3) + x^+(n-1)
x^−n = x^−(n-3) − x^−(n-2)

lua> B, C = {1,1,4}, {0,-2,3}
lua> for n=4,40 do B[n]=B[n-3]+B[n-1]; C[n]=C[n-3]-C[n-2] end
lua> for n=30,36 do print(n, B[n], C[n]) end

Code:

30      95545   406
31      140028  403
32      205221  -870
33      300766  3
34      440794  1273
35      646015  -873
36      946781  -1270

For n=33, linear term coefficient is relatively tiny. Also, (α^n, β^n, γ^n) well separated.

lua> surd(B[33], 33) -- estimated Psi
1.4655712318782408

--

lua> B[33]^2-2*C[33], C[33]^2-2*B[33] --Graeffe Root Squaring
90460186750      -601523

Free42: (90460186750)^(1/66)
1.465571231876768025024430123337594

Free42: (90460186750^2 + 2*601523)^(1/132) // "square" again
1.465571231876768026656731225219939

[Albert Chan]

ψ^3 = (ψ^2 + 1)
ψ^4 = (ψ^2 + 1)*ψ = (ψ^2 + 1) + ψ = (1*ψ^2 + ψ + 1)
ψ^5 = (ψ^2 + 1)*ψ^2 = (ψ^2 + ψ + 1) + ψ^2 = (2*ψ^2 + ψ + 1)
ψ^6 = (ψ^2 + 1)^2 = (ψ^2 + ψ + 1) + 2*ψ^2 + 1 = (3*ψ^2 + ψ + 2)
...
ψ^n = S(n-1)*ψ^2 + S(n-3)*ψ + S(n-2)
ψ^(2n) = S(2n-1)*ψ^2 + S(2n-3)*ψ + S(2n-2)
ψ^(3n) = S(3n-1)*ψ^2 + S(3n-3)*ψ + S(3n-2)

ψ^(3n) = k1*ψ^(2n) + k2*ψ^n + k3

k1 = round(ψ^n), for n ≥ 6      → ψ ≈ k1^(1/n)
k3 = (αβγ)^n = 1

For big n, S(3n-1)/S(3n-3) is better estimate for ψ^2, than S(2n-1)/S(2n-3).
Thus, we expected k2 ≠ 0. But, there are exceptions, when n is small.

n = 1: ψ^3 = ψ^2 + 1      // n=1 --> k2=0

n=11: S(3n-1)/S(3n-3) = 85626/39865 = (67*1278)/(67*595) = 1278/595 = S(2n-1)/S(2n-3)

ψ^33 = 67*ψ^22 + 1      // n=11 --> k2=0

ψ^11 = 67 + 1/ψ^22 = 67 + 1/(67 + 1/ψ^22)^2

(08-09-2022 10:45 AM)Gerson W. Barbosa Wrote:  Or, for more digits,

(67 + 1/(67^2 + 2/(67 + 3/(2×67^2 + 11/(3×67)))))^(1/11)
= 1.46557123187676802665673122475

This also explained why ψ^33 is very close to integer.
With cubing of roots, "linear" term is still small.

(-ψ^11)^3 + 67*(-ψ^11)^2 + 1 = 0

lua> b, c = 67, 0
lua> b*(b*b-3*c)+3, c*(c*c-3*b)+3
300766      3

(-ψ^33)^3 + 300766*(-ψ^33)^2 + 3*(-ψ^33) + 1 = 0

ψ^33 = 300766 - 3/ψ^33 + 1/ψ^66 ≈ 300766 - 3/300766 + 1/300766^2 ≈ 300765.9999900255