How can I do an efficient calculation / approximation of a set of values of a reciprocal function with as few divisions as possible?

Question

I have a function of the form $$f(t) = \frac{a}{bt+c}$$ Which I want to calculate (or numerically approximate) with regular steps on an interval, let us say $$t = t_0 + (t_1-t_0)\frac{k}{N}, \forall k \in\{0,\cdots,N\}$$

The thing is I want to avoid numerical divisions as much as possible as those are slow to do with my computer. Multiplications and additions are okay. But of course I want to get away with a few of those as well (if I can!).

Answer 1

To simplify, let's take $a=b=c=1$ and $t_0=0$ and $t_1=1$.

Define $s_k = \frac{k}N$. You want to evaluate $f$ at $N+1$ points $s_0,\dots,s_{N}$. To beat the direct computation, you need to do it in less than $N+1$ divisions. You don't seem to care so much about the number of additions or multiplications.

$$\begin{split} f(s_{k+1})&=\frac 1 {1+\frac {k+1} N}\\ &=\frac 1 {1+\frac {k} N + \frac{1}N}\\ &=\frac 1 {1+\frac {k} N}\cdot \frac {1}{1+\frac{\frac 1N}{1+\frac{k}{N}}}\\ &=f(s_k)\cdot \frac {1}{1+\frac {f(s_k)}N}\\ &=f(s_k)\cdot \sum_{p\geq 0}(-1)^p\left( \frac {f(s_k)}N\right)^p \,\,\,\,\,\,\,\,\,\text{(1)} \end{split}$$ Truncating the sum above to the $m$-th term yields an error of $$f(s_k)\left( \frac {f(s_k)}N\right)^{m+1}\leq \frac 1 {N^{m+1}}$$ So an algorithm to estimate to a given precision $\varepsilon>0$

• Initialization: Compute $f(s_0)$ [1 division] and $\alpha=\frac 1 N$ [1 division]
• Computation of $f(s_{k+1})$ given $f(s_k)$ and $\alpha$: [no division] $$\left\{ \begin{split} f(s_{k+1})&\simeq f(s_k)\cdot \sum_{p= 0}^m(-1)^p\alpha^p f(s_k)^p\\ \text{where } & \alpha^{m+1}<\varepsilon \end{split} \right.$$

This algorithm will make you do many more additions and multiplications than the direct computation, but you seemed fine with it.

Answer 2

I found a way to do this. We have the function : $$t\to f(t) = \frac{a}{(tb+c)}$$ Consider the derivative of our function $$\frac{\partial f}{\partial t} = \frac{-ab}{(tb+c)^2}$$ Now at first glance this may seem to not help us very much , as we also get a square on top of the reciprocal. There still is a division involved.

But assume that we at some $t_k$ have an approximation to $f(t_k)$.

We realize we can rewrite the derivative as:

$$\frac{\partial f}{\partial t} = -\frac{b}a f(t)^2$$ Now we can now mimic the Euler forward scheme, taking small steps in the direction of the derivative.

$$f(t_k+t_\Delta) \approx f(x_k) + t_\Delta f'(x_k)$$

We see that with the substitution $k=-b/a$ (which can be pre-calculated as well) what is required to calculate by using such a scheme is two multiplications and one addition per new value (after the first).

Here is a sample run for a randomized set of {a,b,c}: