solve_recursive_sequence.sf

#!/usr/bin/ruby

# Find a recursive formula for a given sequence, using the "number wall" method.

# Inspired by the following Mathologer video:
#   Secrets of the lost number walls
#   https://yewtu.be/watch?v=NO1_-qptr6c

func solve_recursive_sequence(seq) {

    var x = Poly(1)

    var a = seq.len.of(1)
    var b = seq.map_cons(2, {|a,b|
        b - a*x
    })

    loop {

        var c = []
        var all_zero = true

        for i in (1 .. b.end-1) {

            b[i]   \\ next
            a[i]   \\ next
            b[i-1] \\ next
            b[i+1] \\ next

            a[i].is_zero && next

            var entry = (b[i-1]*b[i+1] - b[i].sqr)/(-a[i])

            entry.is_nan && next

            if (all_zero) {
                entry.is_zero || (all_zero = false)
            }

            c[i-1] = entry
        }

        if (all_zero) {

            var p = b.first{ defined(_) && .kind_of(Poly) }
            var cf = p.coeffs
            var d  = cf.pop.tail
            var fc = (Poly(cf)/-d -> coeffs)

            var h = Hash(fc.flat...)
            return (^p.degree -> map { h{_} \\ 0 }.flip)
        }

        (a, b) = (b.slice(1,-1), c)
    }
}

func pretty(ker) {
    "a(n) = " + (^ker->flip.grep{ker[_] != 0}.map{|k|
        (ker[k] == 1 ? '' : "#{ker[k].as_rat}*") + "a(n - #{k+1})"
    }.flip.join(' + ') || '0')
}

say pretty(solve_recursive_sequence(20.of{.fib}))                #=> a(n) = a(n-1) + a(n-2)
say pretty(solve_recursive_sequence(20.of{.fib(3)}))             #=> a(n) = a(n-1) + a(n-2) + a(n-3)
say pretty(solve_recursive_sequence(20.of{.fib(4)}))             #=> a(n) = a(n-1) + a(n-2) + a(n-3) + a(n-4)
say pretty(solve_recursive_sequence(20.of{.sqr}))                #=> a(n) = 3*a(n-1) + -3*a(n-2) + a(n-3)
say pretty(solve_recursive_sequence(20.of{.faulhaber(2)}))       #=> a(n) = 4*a(n-1) + -6*a(n-2) + 4*a(n-3) + -1*a(n-4)

assert_eq(
    solve_recursive_sequence([1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, -1, 1, 1, 1, 1, 1, 2, -1, -3, 1, 1, 1, 1, 2, 5, -3, -7, 1, 1, 1, 0, 5, 13, -7, -15, 1, 1, 0, -5, 13, 33, -15, -31, 1, 2, -5, -23, 33, 81, -31, -63, 2, 9, -23, -79, 81, 193, -63, -128, 9, 41, -79, -239, 193, 449, -128, -265, 41, 161, -239]),
    [0,-1,0,-1,0,-1,0,1,0,1,0,1],
)