Variant of the Secant Method

[ttw]

The order of the secant method for finding zeros of a real function can be improved by using spline-like interpolants. This gives a method which has order just under 2 and uses only 1 function evaluation per step.

https://arxiv.org/pdf/2012.04248.pdf

[Valentin Albillo]

(12-09-2020 04:33 AM)ttw Wrote:  The order of the secant method for finding zeros of a real function can be improved by using spline-like interpolants. This gives a method which has order just under 2 and uses only 1 function evaluation per step.

https://arxiv.org/pdf/2012.04248.pdf

Interesting. Perhaps you might consider implementing it yourself for some HP calc and telling us how it compares in terms of speed and memory versus plain old Newton's method {two evaluations of f(x) per iteration), which can be implemented in as few as 34 steps even in the simple HP-10C model (fewer in RPN models with subroutine capability like the HP-11C).

Anyway, thanks for sharing and take care.
V.

[Albert Chan]

Let's try solving x, for sinc(x) = 1/6.5

lua> function f(x) return sin(x)/x - 1/6.5 end
lua> function secant(x1,y1,x2,y2) return x1-y1/(y2-y1)*(x2-x1) end
lua> x1, x2 = 2.711, 2.712
lua> y1, y2 = f(x1) , f(x2)
lua> x3 = secant(x1,y1, x2,y2)
lua> y3 = f(x3)
lua> x3
2.7113129474671775

-- Inverse Interpolation, assumed f(x) is strictly increasing or descreasing around root
-- Fit x = g(y), then evaluate g(0)

lua> x12 = (x2-x1)/(y2-y1)
lua> x23 = (x3-x2)/(y3-y2)
lua> x123 = (x23-x12)/(y3-y1)
lua> x3 + (0-y3)*(x23 + (0-y2)*x123)
2.711312910356982

-- Aikten interpolation, as 1 liner: inner secant = x4, outer secant = x5
-- y1 x1
-- y2 x2 x3
-- y3 x3 x4 x5

lua> secant(secant(x3,y3, x1,y1),y3, x3,y2)
2.711312910356982

-- improved secant, interpolate for slope

lua> y23 = (y3-y2)/(x3-x2)
lua> y12 = (y2-y1)/(x2-x1)
lua> y123 = (y23-y12)/(x3-x1)
lua> slope = y23 + (x3-x2)*y123
lua> x3 - y3/slope
2.7113129103569826

-- Newton's method, actual slope

lua> x, y = x3, y3
lua> slope = (x*cos(x) - sin(x))/x^2
lua> x - y/slope
2.711312910356983

For reference, I used ASINC, for Free-42 Decimal.
Comparing apples to apples, input is float(1/6.5) = 0x1.3b13b13b13b14p-3

5542891849071380 2 55 X^Y ÷ XEQ "ASINC"

2.71131291035698331469084928874754

All 3 methods gives excellent results, doubling accuracy of x3

[Albert Chan]

I extend previous post, with another round of iteration, using all points.
(First iteration matched previous post, done in double precision)

inverse-interpolation:
2.71131291035698220988289207111279
2.71131291035698331469084928874683

improved secant:
2.71131291035698244759230088334763
2.71131291035698331469084928874750

newton's method:
2.71131291035698307703744887554087
2.71131291035698331469084928874753

[Namir]

I implemented the Generalized Secant algorithm in Excel VBA. I used the same order as in the article's example (divided differences for first and second derivatives). The algorithm did better than Newton's method in most cases. The Generalized Secant algorithm needs two initial guesses, while Newton's method needs one. It made a difference in some test if I gave Newton's method the second guess OR averaged the two guesses. The closer of these two values to the neighboring root gave better results.

The algorithm is very interesting!!!

Thanks!

Namir

[Namir]

I was able to implement in Excel VBA a version of the Generalized Secant algorithm using divided differences for the first, second and third derivatives.

Here is the VBA listing for the version that uses two initial guesses and divided differences for the first and second derivatives.

Code:




Cell Mapping
===========
A1 : "X0"
A2 : value for X0
A3 : "X1"
A4 : value for X1
A5 : "Tolerance"
A6 : value for tolerance
A7 : "F(x)"
A8: String containing expression for f(x)=0, such as exp(x)-3*x^2.
B1 : "X"
C1 : "F(X)"  (for Generalized Secant method)
D1 : "X"
E1 : "F(X)"  (for Newton's method, used for comparison)

VBA Listing
===========

Function Fx(ByVal sFx As String, ByVal X As Double) As Double
  sFx = Replace(sFx, "EXP", "!")
  sFx = Replace(sFx, "X", "(" & CStr(X) & ")")
  sFx = Replace(sFx, "!", "EXP")
  Fx = Evaluate(sFx)
End Function

Sub genSecant2()
  Dim sFx As String, X0 As Double, X1 As Double, X2 As Double, X3 As Double
  Dim Toler As Double, Df1 As Double, Df2 As Double
  Dim Row As Integer, h As Double
  Dim Fx0 As Double, Fx1 As Double, Fx2 As Double, Diff As Double
  
  Range("B2:X1000").Clear
  
  X0 = [A2].Value
  X1 = [A4].Value
  Toler = [A6].Value
  sFx = [A8].Value
  sFx = Replace(sFx, " ", "")
  sFx = UCase(sFx)
  
  Fx0 = Fx(sFx, X0)
  Fx1 = Fx(sFx, X1)
  Df1 = Fx0 / (X0 - X1) + Fx1 / (X1 - X0)
  X2 = X1 - Fx1 / Df1
  Fx2 = Fx(sFx, X2)
  Diff = X2 - X1
  
  Row = 2
  Cells(Row, 2) = X2
  Cells(Row, 3) = Fx(sFx, X2)
  Row = Row + 1
  Do Until Abs(Diff) < Toler
    Df1 = Fx2 / (X2 - X1) + Fx1 / (X1 - X2)
    Df2 = Fx2 / (X2 - X1) / (X2 - X0) + Fx1 / (X1 - X2) / (X1 - X0) + Fx0 / (X0 - X2) / (X0 - X1)
    Diff = Fx2 / (Df1 + Df2 * (X2 - X1))
    X3 = X2 - Diff
    X0 = X1
    Fx0 = Fx1
    X1 = X2
    Fx1 = Fx2
    X2 = X3
    Fx2 = Fx(sFx, X3)
    Cells(Row, 2) = X3
    Cells(Row, 3) = Fx2
    Row = Row + 1
  Loop
  
  ' Newton's method
  X1 = [A4].Value
  Fx1 = Fx(sFx, X1)
  Row = 2
  Diff = 2 * Toler
  Do Until Abs(Diff) < Toler
    h = 0.01 * (1 + Abs(X1))
    Diff = h * Fx1 / (Fx(sFx, X1 + h) - Fx1)
    X1 = X1 - Diff
    Fx1 = Fx(sFx, X1)
    Cells(Row, 4) = X1
    Cells(Row, 5) = Fx1
    Row = Row + 1
  Loop
  
End Sub

Output
======

Columns B and C show the refined guesses and their function values for the Generalized Secant Algorithm.
Columns D and E show the refined guesses and their function values for Newton's method.

Here is the version of the VBA code that uses three initial guesses and the divided differences for the first three derivatives:

Code:


Cell Mapping
===========
A1 : "X0"
A2 : value for X0
A3 : "X1"
A4 : value for X1
A5 : "X2"
A6 : value for X2
A7 : "Tolerance"
A8 : value for tolerance
A9 : "F(x)"
A10: String containing expression for f(x)=0
B1 : "X"
C1 : "F(X)"  (for Generalized Secant method)
D1 : "X"
E1 : "F(X)"  (for Newton's method, used for comparison)

VBA Listing
===========

Function Fx(ByVal sFx As String, ByVal X As Double) As Double
  sFx = Replace(sFx, "EXP", "!")
  sFx = Replace(sFx, "X", "(" & CStr(X) & ")")
  sFx = Replace(sFx, "!", "EXP")
  Fx = Evaluate(sFx)
End Function

Sub genSecant3()
  Dim sFx As String, X0 As Double, X1 As Double, X2 As Double, X3 As Double, X4 As Double
  Dim Toler As Double, Df1 As Double, Df2 As Double, Df3 As Double
  Dim Row As Integer, h As Double
  Dim Fx0 As Double, Fx1 As Double, Fx2 As Double, Fx3 As Double, Diff As Double
  
  Range("B2:X1000").Clear
  
  X0 = [A2].Value
  X1 = [A4].Value
  X2 = [A6].Value
  Toler = [A8].Value
  sFx = [A10].Value
  sFx = Replace(sFx, " ", "")
  sFx = UCase(sFx)
  
  Fx0 = Fx(sFx, X0)
  Fx1 = Fx(sFx, X1)
  Fx2 = Fx(sFx, X2)
  Df1 = Fx0 / (X0 - X1) + Fx1 / (X1 - X0)
  Df2 = Fx2 / (X2 - X1) / (X2 - X0) + Fx1 / (X1 - X2) / (X1 - X0) + Fx0 / (X0 - X2) / (X0 - X1)
  X3 = X2 - Fx2 / (Df1 + Df2 * (X2 - X1))
  Fx3 = Fx(sFx, X3)
  Diff = X3 - X2
  
  Row = 2
  Cells(Row, 2) = X3
  Cells(Row, 3) = Fx3
  Row = Row + 1
  Do Until Abs(Diff) < Toler Or Row > 100
    Df1 = Fx3 / (X3 - X2) + Fx2 / (X2 - X3)
    Df2 = Fx3 / (X3 - X1) / (X3 - X2) + Fx2 / (X2 - X1) / (X2 - X3) + Fx1 / (X1 - X2) / (X1 - X3)
    Df3 = Fx3 / (X3 - X2) / (X3 - X1) / (X3 - X0) + _
          Fx2 / (X2 - X3) / (X2 - X1) / (X2 - X0) + _
          Fx1 / (X1 - X3) / (X1 - X2) / (X1 - X0) + _
          Fx0 / (X0 - X3) / (X0 - X2) / (X0 - X1)
    Diff = Fx3 / (Df1 + Df2 * (X3 - X2) + Df3 * (X3 - X2) * (X3 - X1))
    X4 = X3 - Diff
    X0 = X1
    Fx0 = Fx1
    X1 = X2
    Fx1 = Fx2
    X2 = X3
    Fx2 = Fx3
    X3 = X4
    Fx3 = Fx(sFx, X4)
    Cells(Row, 2) = X4
    Cells(Row, 3) = Fx3
    Row = Row + 1
  Loop
  
  ' Newton's method
  X1 = [A4].Value
  Fx1 = Fx(sFx, X1)
  Row = 2
  Diff = 2 * Toler
  Do Until Abs(Diff) < Toler
    h = 0.01 * (1 + Abs(X1))
    Diff = h * Fx1 / (Fx(sFx, X1 + h) - Fx1)
    X1 = X1 - Diff
    Fx1 = Fx(sFx, X1)
    Cells(Row, 4) = X1
    Cells(Row, 5) = Fx1
    Row = Row + 1
  Loop
  
End Sub

Output
======

Columns B and C show the refined guesses and their function values for the Generalized Secant Algorithm.
Columns D and E show the refined guesses and their function values for Newton's method.

Use different worksheets for each of the above methods as the input values in column A are not the same.

I wrote the above code to make it easy to ready. A programmer can optimize the above code by using arrays, matrices, and loops.

[Albert Chan]

Quadratic interpolation is simply 2 linear interpolations, followed by linear interpolation.

x3 y3
x2 y2 y23
x1 y1 y13 y123

Similarly, cubic interpolation is 2 quadratic interpolations, followed by linear interpolation.

x4 y4
x3 y3
x2 y2 y234
x1 y1 y134 y1234

Using this recursive nature of polynomial interpolation, I get this neat Lua code

Code:

function interp(x1,y1,x2,y2,x3,y3,x4,y4)    -- estimate y = p(0)
    if x3 then
        y2 = interp(x2,y2,x3,y3,x4,y4)
        y1 = interp(x1,y1,x3,y3,x4,y4)
    end
    local y, p = y2, x2/(x2-x1)
    if p > 0.5 then y, p = y1, x1/(x2-x1) end
    return y - p*(y2-y1)
end

Extending interp(...) is trivial, by adding (x5,y5), ...
The down-side is each bump up in polynomial degree, doubled the work.
For high degree fit, divided-difference style is better, because it can cache previous results.

OTTH, assuming we have convergence, old points are very, very far away.
Fitting those probably may not help much ...

---

For slope of improved secant, we can also consider this an interpolation problem.
Example, for a cubic fit, we generate divided-difference table:

x4 y4
x3 y3 y34
x2 y2 y24 y234
x1 y1 y14 y123 y1234

p(x) = y4 + (x-x4)*(y34 + (x-x3)*(y234 + (x-x2)*y1234)) = y4 + (x-x4)*M

p'(x) = M + (x-x4)*M'       → p'(x4) = M

Basing all values relative to (x4, y4). Example X3 = x3-x4, Y3 = y3-y4 ...
M is just interpolation of individual slopes.

X3      Y3/X3
X2      Y2/X2
X1      Y1/X1      M

Again, extending this is simply adding more slopes to fit

Code:

function secant(x1,y1,x2,y2,x3,y3,x4,y4)    -- estimate root of p(x)
    if not x3 then return interp(y1,x1,y2,x2) end
    if x4 then x3,y3,x4,y4 = x4,y4,x3-x4,(y3-y4)/(x3-x4) end
    x1, y1, x2, y2 = x1-x3, y1-y3, x2-x3, y2-y3
    return x3 - y3 / interp(x1,y1/x1, x2,y2/x2, x4,y4)    
end

Example:

lua> function f(x) return sin(x)/x - 1/6.5 end
lua> x1, x2 = 2.711, 2.712
lua> y1, y2 = f(x1), f(x2)
lua> x3 = secant(x1,y1, x2,y2)
lua> y3 = f(x3)

lua> interp(y1,x1, y2,x2, y3,x3) -- inverse interpolate, x = g(0)
2.711312910356982
lua> secant(x1,y1, x2,y2, x3,y3) -- improved secant, x = x3 - y3/slope
2.7113129103569826

[Albert Chan]

Let's compare inverse-interpolation vs improved-secant (code from previous post)

Code:

function test(f, x0, x1, n)   -- n = 1 (linear) to 3 (cubic)
    local y0, y1 = f(x0), f(x1)
    local t0 = {y0,x0, y1,x1} -- inverse interpolate
    local t1 = {x0,y0, x1,y1} -- improved secant
    n = (n or 3)*2 + 1
    for i = 4,18,2 do         -- i = table length
        x0 = interp(unpack(t0, math.max(1,i-n)))
        x1 = secant(unpack(t1, math.max(1,i-n)))
        t0[i+1], t0[i+2] = f(x0), x0
        t1[i+1], t1[i+2] = x1, f(x1)
        print(i/2, x0, x1)
    end
end

lua> function f(x) return x^3 - 8 end

lua> test(f, 5, 4, 2) -- https://arxiv.org/pdf/2012.04248.pdf, Table 2
2      3.081967213114754        3.081967213114754
3      2.3945608933295457      2.2862188297178117
4      2.0980163342039635      2.010344209437878
5      2.0088865345029276      1.99979593345267
6      2.0001129292461193      2.0000000722313933
7      2.0000000388749792      2.000000000000015
8      2.0000000000000164      2
9      2                                    2

Last column is improved secant, fitting data up to quadratic (slope by linear interpolation)
Middle column is inverse-interpolation for x, at y = 0, also fitting up to quadratic.

Improved secant uses the last point as the "base", then does correction.
Inverse-interpolation is not as good, probably due to putting equal weight to its points.

Let's try again, but with up to cubic fit.

lua> test(f, 5, 4, 3)
2      3.081967213114754        3.081967213114754
3      2.3945608933295457      2.2862188297178117
4      2.082744738020821        2.034337291023909
5      2.0058425366743506      2.000576313421517
6      2.0000383915262443      2.000000166004798
7      2.0000000023200664      2.0000000000000138
8      1.9999999999999998      2
9      2                                    2

Based on above, anything more than a cubic fit probably not worth the trouble.
For example, lets see how Newton's method does (actual slope)

lua> function df(x) return 3*x*x end
lua> x = 4
lua> for i=2, 9 do x = x - f(x)/df(x); print(i, x) end
2      2.833333333333333
3      2.221068819684737
4      2.0212735368091126
5      2.0002231146078984
6      2.0000000248863623
7      2.0000000000000004
8      2
9      2

[Namir]

I went back to my Excel VBA code and added statements to implement Halley's method. The Generalized Secant method did better than Newton's method but ws behind Halley's method in the number of steps. Of course results depend on the function tested and the initial guess.

Namir

[ttw]

Methods with high convergence can have problems numerically; they also tend to have a smaller domain of convergence. Multi-point methods seem to be a bit more stable than high order methods.

[Albert Chan]

(12-12-2020 02:41 PM)ttw Wrote:  Methods with high convergence can have problems numerically; they also tend to have a smaller domain of convergence.

This might not always true.
In fact, recently we come across opposite case, LambertW function implementation.

Solving f(x) = x*e^x - a, with Newton's method is very unstable.

Higher convergence method, say Halley, it becomes "safe"

P.S. the next post, I had suggested another way, to suppress the exponential.
This version, Newton's method is fast, and safe.

>>> from mpmath import *
>>> a = -1e99
>>> findroot(lambda x: x + ln(x/a), 200+3j) # bad guess
mpc(real='222.55067066150301', imag='3.1275404175517207')
>>> lambertw(a)
mpc(real='222.55067066150301', imag='3.127540417551721')

[Albert Chan]

(12-11-2020 04:34 PM)Albert Chan Wrote:  Quadratic interpolation is simply 2 linear interpolations, followed by linear interpolation.

x3 y3
x2 y2 y23
x1 y1 y13 y123

It may seems we need 3 secant root steps for 3-points fit, in a loop, we only need 2

x2 y2
x3 y3 y23
x1 y1 y12 y123

Next iteration y23 → y12, can be reused.

Secant setup code

Code:

function newx(a,fa, b,fb) return b - fb/(fb-fa)*(b-a) end

function secant_init(f,a,b,eps,verbal)
    eps = eps or 1e-9
    local f2 = verbal and function(x) print(x); return f(x) end or f
    local fa = f2(a); if fa==0 then return nil, a end
    local fb = f2(b); if fb==0 then return nil, b end
    if signbit(eps) and abs(fa)<abs(fb) then
        a,fa, b,fb = b,fb, a,fa
        if verbal then print('-- points swapped') end
    end
    return f2, a,fa, b,fb, abs(eps)
end

Code:

function S.secant2(f,a,b,eps,verbal,c)  -- 3 pt fit
    local fa,fb,fc,ab,bc
    f, a,fa, b,fb, eps = secant_init(f,a,b,eps,verbal)
    if not f then return a,a end
    ab = newx(a,fa, b,fb)
    c = c or ab
    fc = f(c); if fc==0 then return c,c end
    while abs(b-c) > eps do
        bc, ab = ab, newx(b,fb, c,fc)
        a,fa, b,fb, c = b,fb, c,fc, newx(bc,fa, ab,fc)
        fc = f(c); if fc==0 then return c,c end
    end
    return finite(newx(b,fb, c,fc)) or (b+c)/2
end

lua> f = fn'x: x^3 - 8'
lua> x = S.secant2(f, 5, 4, 1e-9, true)
5
4
3.081967213114754
2.3945608933295457
2.0980163342039635
2.0088865345029276
2.0001129292461193
2.0000000388749792
2.0000000000000164
2

Result matched post #8

---

For completeness, here is S.secant3 (OP "improved" secant)

Code:

function S.secant3(f,a,b,eps,verbal,c)  -- slope extrapolated
    local fa,fb,fc
    f, a,fa, b,fb, eps = secant_init(f,a,b,eps,verbal)
    if not f then return a,a end
    c = c or newx(a,fa, b,fb)
    fc = f(c); if fc==0 then return c,c end
    while abs(b-c) > eps do
        local ca,cb = c-a,c-b
        local slope = newx((fc-fa)/ca,ca, (fc-fb)/cb,cb)
        a,fa, b,fb, c = b,fb, c,fc, c-fc/slope
        fc = f(c); if fc==0 then return c,c end
    end
    return finite(newx(b,fb, c,fc)) or (b+c)/2
end

lua> x = S.secant3(f, 5, 4, 1e-9, true)
5
4
3.081967213114754
2.2862188297178117
2.010344209437878
1.99979593345267
2.0000000722313933
2.000000000000015
2

[Thomas Klemm]

These functions implement the described secant method in Python for orders \(k \in \{1, 2, 3\}\):

Code:

def f(x):
    print(x)
    return x**3 - 8


def secant_1(x_4, x_5, ε):
    f_4 = f(x_4)
    f_5 = f(x_5)
    f_45 = (f_4 - f_5) / (x_4 - x_5)

    while ε < abs(f_5):
        x_6 = x_5 - f_5 / f_45
        f_6 = f(x_6)
        f_56 = (f_5 - f_6) / (x_5 - x_6)
        x_5, f_5, f_45 = x_6, f_6, f_56

    return x_5


def secant_2(x_4, x_5, ε):
    f_4 = f(x_4)
    f_5 = f(x_5)
    f_45 = (f_4 - f_5) / (x_4 - x_5)

    x_6 = x_5 - f_5 / f_45
    f_6 = f(x_6)
    f_56 = (f_5 - f_6) / (x_5 - x_6)
    f_456 = (f_45 - f_56) / (x_4 - x_6)

    while ε < abs(f_6):
        x_7 = x_6 - f_6 / (f_56 + (x_6 - x_5) * f_456)
        f_7 = f(x_7)
        f_67 = (f_6 - f_7) / (x_6 - x_7)
        f_567 = (f_56 - f_67) / (x_5 - x_7)
        x_5, x_6 = x_6, x_7
        f_6, f_56, f_456 = f_7, f_67, f_567

    return x_6


def secant_3(x_4, x_5, ε):
    f_4 = f(x_4)
    f_5 = f(x_5)
    f_45 = (f_4 - f_5) / (x_4 - x_5)

    x_6 = x_5 - f_5 / f_45
    f_6 = f(x_6)
    f_56 = (f_5 - f_6) / (x_5 - x_6)
    f_456 = (f_45 - f_56) / (x_4 - x_6)

    x_7 = x_6 - f_6 / (f_56 + (x_6 - x_5) * f_456)
    f_7 = f(x_7)
    f_67 = (f_6 - f_7) / (x_6 - x_7)
    f_567 = (f_56 - f_67) / (x_5 - x_7)
    f_4567 = (f_456 - f_567) / (x_4 - x_7)

    while ε < abs(f_7):
        x_8 = x_7 - f_7 / (f_67 + (x_7 - x_6) * (f_567 + f_4567 * (x_7 - x_5)))
        f_8 = f(x_8)
        f_78 = (f_7 - f_8) / (x_7 - x_8)
        f_678 = (f_67 - f_78) / (x_6 - x_8)
        f_5678 = (f_567 - f_678) / (x_5 - x_8)
        x_5, x_6, x_7 = x_6, x_7, x_8
        f_7, f_67, f_567, f_4567 = f_8, f_78, f_678, f_5678

    return x_7

We can run them using the Decimal library:

Code:

from decimal import Decimal, getcontext
getcontext().prec = 34

Code:

secant_1(Decimal(5), Decimal(4), 1e-32)

5
4
3.081967213114754098360655737704918
2.519552120040923041946117994770945
2.180972989759050190092856538085551
2.037953100909517790045306042203723
2.003198489980016115117431120840932
2.000059872823468592338193344700047
2.000000095647401657566635004047756
2.000000000002863282761488899959229
2.000000000000000000136932773807772
2.000000000000000000000000000000196
2.000000000000000000000000000000000

Code:

secant_2(Decimal(5), Decimal(4), 1e-32)

5
4
3.081967213114754098360655737704918
2.286218829717811307322668037730626
2.010344209437878312641529731720143
1.999795933452669925783583536567984
2.000000072231393330599606713662298
2.000000000000015319238844912588532
2.000000000000000000000000018934481
2.000000000000000000000000000000000

Code:

secant_3(Decimal(5), Decimal(4), 1e-32)

5
4
3.081967213114754098360655737704918
2.286218829717811307322668037730626
2.034337291023909027923796138229571
2.000576313421516741692818211990178
2.000000166004797850203846960058336
2.000000000000013778794929746133383
2.000000000000000000000000000094928
2.000000000000000000000000000000000


These are the corresponding programs for the HP-42S:

Code:

00 { 76-Byte Prgm }
01▸LBL "SEC-1"
02 STO 02 # ε
03 R↓
04 STO 01 # x_5
05 R↓
06 STO 00 # x_4
07 XEQ F
08 STO 03 # f_4
09 RCL 01 # x_5
10 XEQ F
11 STO 04 # f_5
12 RCL 03
13 RCL- 04
14 RCL 00
15 RCL- 01
16 ÷
17 STO 05 # f_45 = (f_4 - f_5) ÷ (x_4 - x_5)
18▸LBL 00
19 RCL 02
20 RCL 04
21 ABS
22 X<Y?   # abs(f_5) < ε
23 GTO 01
24 RCL 01
25 RCL 04
26 RCL 05
27 ÷
28 -
29 STO 06 # x_6 = x_5 - f_5 ÷ f_45
30 XEQ F
31 X<> 04 # f_6 <> f_5
32 RCL- 04
33 RCL 06
34 X<> 01 # x_6 <> x_5
35 RCL- 01
36 ÷
37 STO 05 # f_45 = f_56 = (f_5 - f_6) ÷ (x_5 - x_6)
38 GTO 00
39▸LBL 01
40 RCL 01 # x_5
41 RTN
42▸LBL F
43 STOP
44 3
45 Y↑X
46 8
47 -
48 END

Code:

00 { 127-Byte Prgm }
01▸LBL "SEC-2"
02 STO 02 # ε
03 R↓
04 STO 01 # x_5
05 R↓
06 STO 00 # x_4
07 XEQ F
08 STO 03 # f_4
09 RCL 01 # x_5
10 XEQ F
11 STO 04 # f_5
12 RCL 03
13 RCL- 04
14 RCL 00
15 RCL- 01
16 ÷
17 STO 05 # f_45 = (f_4 - f_5) ÷ (x_4 - x_5)
18 RCL 01
19 RCL 04
20 RCL 05
21 ÷
22 -
23 STO 06 # x_6 = x_5 - f_5 ÷ f_45
24 XEQ F
25 STO 07 # f_6
26 RCL 04
27 RCL- 07
28 RCL 01
29 RCL- 06
30 ÷
31 STO 08 # f_56 = (f_5 - f_6) ÷ (x_5 - x_6)
32 RCL 05
33 RCL- 08
34 RCL 00
35 RCL- 06
36 ÷
37 STO 09 # f_456 = (f_45 - f_56) ÷ (x_4 - x_6)
38▸LBL 00
39 RCL 02
40 RCL 07
41 ABS
42 X<Y?   # abs(f_6) < ε
43 GTO 01
44 RCL 06
45 RCL 07
46 RCL 06
47 RCL- 01
48 RCL× 09
49 RCL+ 08
50 ÷
51 -
52 STO 10 # x_7 = x_6 - f_6 ÷ ((x_6 - x_5) × f_456 + f_56)
53 XEQ F
54 X<> 07 # f_7 <> f_6
55 RCL- 07
56 RCL 06
57 RCL- 10
58 ÷
59 X<> 08 # f_56 <> f_67 = (f_6 - f_7) ÷ (x_6 - x_7)
60 RCL- 08
61 RCL 01
62 RCL- 10
63 ÷
64 STO 09 # f_456 = f_567 = (f_56 - f_67) ÷ (x_5 - x_7)
65 RCL 10
66 X<> 06
67 STO 01 # x_5, x_6 = x_6, x_7
68 GTO 00
69▸LBL 01
70 RCL 06 # x_6
71 RTN
72▸LBL F
73 STOP
74 3
75 Y↑X
76 8
77 -
78 END

Code:

00 { 195-Byte Prgm }
01▸LBL "SEC-3"
02 STO 02 # ε
03 R↓
04 STO 01 # x_5
05 R↓
06 STO 00 # x_4
07 XEQ F
08 STO 03 # f_4
09 RCL 01 # x_5
10 XEQ F
11 STO 04 # f_5
12 RCL 03
13 RCL- 04
14 RCL 00
15 RCL- 01
16 ÷
17 STO 05 # f_45 = (f_4 - f_5) ÷ (x_4 - x_5)
18 RCL 01
19 RCL 04
20 RCL 05
21 ÷
22 -
23 STO 06 # x_6 = x_5 - f_5 ÷ f_45
24 XEQ F
25 STO 07 # f_6
26 RCL 04
27 RCL- 07
28 RCL 01
29 RCL- 06
30 ÷
31 STO 08 # f_56 = (f_5 - f_6) ÷ (x_5 - x_6)
32 RCL 05
33 RCL- 08
34 RCL 00
35 RCL- 06
36 ÷
37 STO 09 # f_456 = (f_45 - f_56) ÷ (x_4 - x_6)
38 RCL 06
39 RCL 07
40 RCL 06
41 RCL- 01
42 RCL× 09
43 RCL+ 08
44 ÷
45 -
46 STO 10 # x_7 = x_6 - f_6 ÷ ((x_6 - x_5) × f_456 + f_56)
47 XEQ F
48 STO 11 # f_7
49 RCL 07
50 RCL- 11
51 RCL 06
52 RCL- 10
53 ÷
54 STO 12 # f_67 = (f_6 - f_7) ÷ (x_6 - x_7)
55 RCL 08
56 RCL- 12
57 RCL 01
58 RCL- 10
59 ÷
60 STO 13 # f_567 = (f_56 - f_67) ÷ (x_5 - x_7)
61 RCL 09
62 RCL- 13
63 RCL 00
64 RCL- 10
65 ÷
66 STO 14 # f_4567 = (f_456 - f_567) ÷ (x_4 - x_7)
67▸LBL 00
68 RCL 02
69 RCL 11
70 X<Y?   # abs(f_7) < ε
71 GTO 01
72 RCL 10
73 RCL 11
74 RCL 10
75 RCL- 01
76 RCL× 14
77 RCL+ 13
78 RCL 10
79 RCL- 06
80 ×
81 RCL+ 12
82 ÷
83 -
84 STO 15 # x_8 = x_7 - f_7 ÷ (((x_7 - x_5) × f_4567 + f_567) × (x_7 - x_6) + f_67)
85 XEQ F
86 X<> 11 # f_8 <> f_7
87 RCL- 11
88 RCL 10
89 RCL- 15
90 ÷
91 X<> 12 # f_67 <> f_78 = (f_7 - f_8) ÷ (x_7 - x_8)
92 RCL- 12
93 RCL 06
94 RCL- 15
95 ÷
96 X<> 13 # f_567 <> f_678 = (f_67 - f_78) ÷ (x_6 - x_8)
97 RCL- 13
98 RCL 01
99 RCL- 15
100 ÷
101 STO 14 # f_4567 = f_5678 = (f_567 - f_678) ÷ (x_5 - x_8)
102 RCL 15
103 X<> 10
104 X<> 06
105 STO 01 # x_5, x_6, x_7 = x_6, x_7, x_8
106 GTO 00
107▸LBL 01
108 RCL 10 # x_7
109 RTN
110▸LBL F
111 STOP
112 3
113 Y↑X
114 8
115 -
116 END

From the paper: \(s_1 \doteq 1.618, \; s_2 \doteq 1.839, \; s_3 \doteq 1.928\)
It appears that the benefit of order 3 compared to order 2 is not that big.

Registers

Interestingly the 3 variants share registers and also some of the code.

SEC-1:
00: x_4
01: x_5
02: ε
03: f_4
04: f_5
05: f_45 = (f_4 - f_5) ÷ (x_4 - x_5)
06: x_6 = x_5 - f_5 ÷ f_45

SEC-2:
00: x_4
01: x_5
02: ε
03: f_4
04: f_5
05: f_45 = (f_4 - f_5) ÷ (x_4 - x_5)
06: x_6 = x_5 - f_5 ÷ f_45
07: f_6
08: f_56 = (f_5 - f_6) ÷ (x_5 - x_6)
09: f_456 = (f_45 - f_56) ÷ (x_4 - x_6)
10: x_7 = x_6 - f_6 ÷ ((x_6 - x_5) × f_456 + f_56)

SEC-3:
00: x_4
01: x_5
02: ε
03: f_4
04: f_5
05: f_45 = (f_4 - f_5) ÷ (x_4 - x_5)
06: x_6 = x_5 - f_5 ÷ f_45
07: f_6
08: f_56 = (f_5 - f_6) ÷ (x_5 - x_6)
09: f_456 = (f_45 - f_56) ÷ (x_4 - x_6)
10: x_7 = x_6 - f_6 ÷ ((x_6 - x_5) × f_456 + f_56)
11: f_7
12: f_67 = (f_6 - f_7) ÷ (x_6 - x_7)
13: f_567 = (f_56 - f_67) ÷ (x_5 - x_7)
14: f_4567 = (f_456 - f_567) ÷ (x_4 - x_7)
15: x_8 = x_7 - f_7 ÷ (((x_7 - x_5) × f_4567 + f_567) × (x_7 - x_6) + f_67)
16: f_8

Example

This is the example from the paper in table 2.
Keep hitting the R/S key to get the next result.

5 ENTER
4 ENTER
1E-32
XEQ "SEC-2"


5
4
3.081967213114754098360655737704918
2.286218829717811307322668037730626
2.010344209437878312641529731720143
1.999795933452669925783583536567984
2.000000072231393330599606713662298
2.000000000000015319238844912588532
2.000000000000000000000000018934481
2.000000000000000000000000000000000