sidef-scripts/Math/function_inverse_mediant_inequality.sf

#!/usr/bin/ruby

# Daniel "Trizen" Șuteu
# Date: 28 July 2018
# https://github.com/trizen

# Compute the inverse of any function by narrowing the search using the mediant inequality:
#   a/b < (a+c)/(b+d) < c/d

# This algorithm, on average, takes more steps to converge than the binary search algorithm.

# See also:
#   https://en.wikipedia.org/wiki/Farey_sequence
#   https://en.wikipedia.org/wiki/Stern%E2%80%93Brocot_tree
#   https://en.wikipedia.org/wiki/Mediant_(mathematics)

func mediant_inverse (n, f, min=0, max=n, prec=192) {

    local Number!PREC = prec.numify

    (min, max) = (max, min) if (min > max)

    var (ln, ld) = min.nude
    var (rn, rd) = max.nude

    loop {
        var m = (ln+rn)/(ld+rd)
        var c = (f(m) <~> n)

        if (c < 0) {
            (ln, ld) = m.nude
        }
        elsif (c > 0) {
            (rn, rd) = m.nude
        }
        else {
            return m
        }
    }
}

say mediant_inverse(  2, {|x| x.exp })     # solution to x for: exp(x) =   2
say mediant_inverse( 43, {|x| x**2  })     # solution to x for:    x^2 =  43
say mediant_inverse(-43, {|x| x**3  })     # solution to x for:    x^3 = -43

# Find the value of x such that Li(x) = 100
say mediant_inverse(100, {|x| Li(x) }, min: 1, max: 1e6, prec: 200)    #=> 488.87190985280753190605086392033334827378018556409