Problem:  You are given an invertable function F of a single variable over a finite interval. This function might take a long time to compute. Write a program that computes a good approximation to the function’s inverse at any point, and that is always fast. You are allowed one arbitrarily long setup computation.

In other words, write a program that takes a given y in the function’s range and quickly computes x so that F(x) = y to a good approximation. And this even though F might take a long time to compute.

The solution is to use function interpolation on the inverse, using some algorithm to compute the inversion table. Any of the usual root finding algorithms (bisection, secant, Brent’s, or – if the derivative is known – Newton’s) could be used for the latter. Since this computation needn’t be very fast we use the bisection method, which is always reliable. (If you use one of the other root finding algorithms, note that F is strictly monotonic because it is invertable, which might simplify the computation.)

Using the bisection algorithm (or any other) by itself will not solve the original problem because it takes even longer to compute than the function. The method shown below, after its long initial computation, computes the inverse at once for any given y.

Before we get to the method first a description of the two algorithms it depends on.


Bisection algorithm for solving F(x) = 0

The accuracy is determined by the constant tol (tolerance). No upper bound has been placed on the number of iterations. In order to avoid an infinite loop tol must not be smaller than the machine precision used for x.

Two initial numbers must be given that bracket the solution
a  ≤  x  ≤  b

We assume the function is invertable over the given domain which means there is only one root.

The program below and the ones given later can be modified for Visual Basic, FORTRAN and any other BASIC-like compiler.

Function F(x As Single) As Single
 Function = x*x - 2   'simple example, the solution in the domain  1 < x < 2  is  Sqr(2)
End Function
'----------------------------

Function Main
 Local i As Long
 Local x, x1, x2 As Single
 Local tol As Single

 x1 = 1               'bracket solution
 x2 = 2
 tol = .000001        'fixed
 i = 0                'count loops, i not used otherwise
 
 Do
  Incr i
  x = (x1 + x2) * .5
  If F(x) = 0 Or Abs(x1 - x2) < tol Then Exit Do
  If Sgn(F(x)) = Sgn(F(x2)) Then
   x2 = x
  Else
   x1 = x
  End If
 Loop
 
 Print
 Print tol; "tolerance"
 Print i; "iterations"
 Print " Solution: "; x
 Print " Actual solution: "; Sqr(2)
 WaitKey$
End Function

Given evaluations of G(y) at evenly spaced y over an interval, find a good approximation of G(y) for any y in the interval by linear interpolation.

This is useful when G takes a lot of time to compute and great accuracy is unnecessary. You compute G at a finite number of values once and for all, then use this table to approximate G at any value any number of times.

The greater the finite number of these values (n) or the shorter the interval over which they are spread (interval), the more accurate the interpolation.

In order to minimize truncation error we put up with a trivial time penalty by using
      ... * n / interval
instead of pre-computing the fraction.

Usually y is considered a function of x, but here the relation is the reverse, because we have in mind an application to the inverse of a function.

In this example the function is G(y) = Sqr(y).  You can change this.

The interpolation will be exact at the y() and less accurate in between. For example, with the fixed parameters in the program as is, where
      n = 50,  y0 = 0,  yn = 5,  G(y) = Sqr(y)
consider  y = 4:
      Ginterpolate(4)  =  2
      G(4)  =  2
With 4 it’s exact, but at  y = 3.95
      Ginterpolate(3.95)  =  1.987340
      G(3.95)  =  1.987461

Global y0, yn As Single            'first and last points
Global interval As Single          'length of interval from y0 to yn
Global n As Long                   'number of divisions plus one
Global x(), y() As Single          'x = G(y),  y() equally spaced from y0 to yn
'----------------------------

' The larger n is (see Sub SetupG) the more accurate the interpolation.

Function G(y As Single) As Single  'defined from y0 to yn
 Function = Sqr(y)                 'simple example
End Function
'----------------------------

Sub SetupG
 n = 50             'fixed

 Dim x(n), y(n) As Single
 Local i As Long

' chose the first and last points so the line through them on the graph 
' does not intersect a zero of the function
 y0 = 0             'first point, fixed
 yn = 10            'last point, fixed
 interval = yn - y0

 For i = 0 To n - 1
  y(i) = y0 + interval * i / n      ' y0 + s*i  where  s = (yn - y0) / n
  x(i) = G(y(i))
 Next
 y(n) = yn          'possibly more accurate than putting in loop
 x(n) = G(yn)
End Sub
'----------------------------

Function Ginterpolate(y As Single) As Single
 Local k As Long
 k = Int((y - y0)*n / interval)     ' Int((y - y0) / s)  where  s = (yn - y0) / n
 Function = x(k) + (y - y(k)) * ( x(k+1) - x(k) ) * n / interval    ' ...*n/interval = .../s
End Function
'----------------------------

Function Main
 Local y As Single
 Local a As String

 SetupG

 Do
  Line Input "Enter a number:  "; a  : If Trim$(a) = "" Then End
  y = Val(a)
  If y < y0 Or y > yn Then    'check range
   Print "Number must be between" + Str$(y0) + " and" + Str$(yn)
   Iterate
  End If
  Print "Interpolated: "; Ginterpolate(y); "     From function: "; G(y)
  Print
 Loop
End Function

For the purpose of demonstrating the idea described in the introduction the particular function to invert below is fast, simply x squared, and has a known simple solution. You can change the function in Function F, then either change or delete the “Actual inverse” in the output.

Global y0, yn As Single            'first and last points
Global interval As Single          'length of interval from y0 to yn
Global n As Long                   'number of interpolation points plus one
Global x(), y() As Single          'x = G(y),  y() equally spaced from y0 to yn
Global tol As Single
'----------------------------

Function F(x As Single) As Single
 Function = x*x    'simple example for the purpose of demonstration
End Function
'----------------------------

Function G(y As Single) As Single  'the inverse of F
 Local x, x1, x2 As Single
 Local y1, y2 As Single
 x1 = y0
 x2 = yn
 Do
  x = (x1 + x2) * .5
  If F(x) = y Or Abs(x1 - x2) < tol Then Exit Do
  If Sgn(F(x) - y) = Sgn(F(x2) - y) Then
   x2 = x
  Else
   x1 = x
  End If
 Loop
 Function = x
End Function
'----------------------------

Sub SetupG               'called only once
 Dim x(n), y(n) As Global Single
 Local i As Long
 interval = yn - y0
 For i = 0 To n - 1
  y(i) = y0 + interval * i / n
  x(i) = G(y(i))
 Next
 y(n) = yn
 x(n) = G(yn)
End Sub
'----------------------------

Function Ginterpolate(y As Single) As Single
 Local k As Long
 k = Int((y - y0)*n / interval)
 Function = x(k) + (y - y(k)) * ( x(k+1) - x(k) ) * n / interval
End Function
'----------------------------

Function Main
 Local y As Single
 Local a As String
 y0 = 1                                      '<--- change to suit
 yn = 10                                     '<--- change to suit
 tol = .000001                               'must be no smaller than precision of x
 SetupG
 Print " Range:  "; Format$(y0);" to ";Format$(yn)
 Print " Number of interpolation points:  "; Format$(n)
 Print " Tolerance:  ";  Format$(tol,"0.######")
 Print
 Do
  Line Input " Enter a number:  "; a  : If Trim$(a) = "" Then End
  y = Val(a)
  If y < y0 Or y > yn Then                   'check range
   Print " Number must be between" + Str$(y0) + " and" + Str$(yn)
   Iterate
  End If
  Print " Interpolate:"; Ginterpolate(y)
  Print " Our inverse:"; G(y)                '<--- not used in practice
  Print " Actual inverse:"; Sqr(y)           '<--- not used in practice
  Print
 Loop
End Function