123456789101112131415161718192021222324252627282930313233343536373839404142434445464748 |
- #!/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
|