[Math] Scaling numbers cleverly to prevent arithmetic overflows or rounding errors

arithmetic

Let's say I have an integer position $p$ inside some integer range $r$, such that

$$
0 \le p < r \;\;\;\textrm{with}\;p, r \in \mathbb{N}
$$

I need to scale a position $p$ from a small range where $0<r<2^{31}$ to the corresponding position $P$ in a big range $0<R<2^{63}$ and back (in programming terms: I have one scenario where my position and range are 32-bit integers and another where my position and range are 64-bit integers).

So, essentially all I need to do is:

$$
P = \lfloor \, p * R / r\rfloor\;\;\;\;\;\textrm{and}\;\;\;\;\
p = \lfloor P * r / R \rfloor
$$
But, for this transformation, I am only allowed to perform 64-bit integer operations.

To scale the small range up to the big range without generating overflows or rounding errors is relatively easy:

$$
P = \lfloor \, p * R / r\rfloor = \lfloor R / r\rfloor * p + \lfloor((R\bmod{r} ) * p) / r\rfloor
$$
(I used $\lfloor\cdot\rfloor$ to indicate the implicit rounding of the integer division operations)

Aside –
a quick explanation of the reasoning behind this equation: Since $p<r$, the first term in the sum will never exceed $2^{63}$. But, doing the division "first" to make this happen, leads to a rounding error. This error is captured by the modulus operation in the second term, which, since $p<r<2^{31}$, can also not exceed $2^{63}$. Thus, this formula allows me to calculate $P$ such that it is correctly rounded down to the nearest integer, only using integer arithmetic with numbers $<2^{63}$.

What I'm having some difficulty with is finding the reverse formula, i.e. for $p = \lfloor P * r / R \rfloor$:

I'm looking for a function $f(P, r, R)$ that calculates $\lfloor P * r / R \rfloor$ while only using 64-bit integer arithmetic, i.e. every division implies an immediate rounding down and all intermediate values must fall within the range $[{-\,2^{63}, +\,2^{63}}[$.

UPDATE:

I've been playing with a recursive solution by partially solving the reverse equation for $p$:

$$ P = \lfloor R/r \rfloor * p + \lfloor ((R\bmod{r})*p)/r \rfloor $$

$$ \lfloor R/r \rfloor * p_{n+1} = P – \lfloor ((R\bmod{r})*p_n)/r \rfloor $$

$$ p_{n+1} = \bigg\lfloor \frac{P}{\lfloor R/r \rfloor} \bigg\rfloor – \bigg\lfloor \frac{\lfloor ((R\bmod{r})*p_n)/r \rfloor} {\lfloor R/r \rfloor} \bigg\rfloor $$

For most values of $r$, $P$, and $R$ it very quickly gives a value that is within $\pm 1$ of the actual answer. But sometimes, it ends up oscillating between two consecutive values and if $R\bmod{r}$ becomes much larger than $R/r$, convergence is either really slow, or sometimes doesn't happen at all.

Can this be improved somehow to guarantee quick convergence to the actual answer? Or maybe it's possible to make this non-recursive?

Best Answer

Disclaimer: This isn't a 100% complete answer, but for the case that I need it for, it works. Maybe someone with additional ideas can improve it...

One can make some progress by looking at different possible special cases:

If $P*r<2^{63}$, there is no problem with overflows to begin with, so we can just calculate the answer directly as $p=\lfloor P*r/R \rfloor$. The only thing we need to do is to perform the check as $P<\lfloor2^{63}/r \rfloor$ instead to allow for a false result without overflow.

To cover the other case, one can do the following transformation:

$$ p=\lfloor P*r/R \rfloor= \bigg\lfloor\frac{P}{R/r}\bigg\rfloor=\bigg\lfloor\frac{P}{\lfloor R/r \rfloor + (R\bmod{r})/r}\bigg\rfloor $$

$$ =\bigg\lfloor\frac{P}{\lfloor R/r \rfloor}-\frac{P(R\bmod{r})/r}{\lfloor R/r \rfloor(R/r)}\bigg\rfloor=\Bigg\lfloor\bigg\lfloor\frac{P}{\lfloor R/r\rfloor}\bigg\rfloor + \frac{P\bmod{\lfloor R/r\rfloor}}{\lfloor R/r\rfloor}-\frac{P(R\bmod{r})}{R\lfloor R/r \rfloor}\Bigg\rfloor $$

$$ =\bigg\lfloor\frac{P}{\lfloor R/r\rfloor}\bigg\rfloor + \Bigg\lfloor\frac{P\bmod{\lfloor R/r\rfloor}}{\lfloor R/r\rfloor}-\frac{P}{R}\frac{R\bmod{r}}{\lfloor R/r \rfloor}\Bigg\rfloor $$

The first term in the sum can be calculated directly without overflow issues. The second half after the '$+$' - sign is a little more difficult. But, looking at it, one realizes the following: The fraction to the left of the '$-$' - sign is always between $\ge 0$ and $<1$. Also: $0 \le P/R < 1$. So, we can look at two different special cases here:

If $R\ge r^2$, the entire part to the right of the '$+$' - sign can only be either $0$ or $-1$.
Otherwise, things could get more complicated...

Luckily, in my application, I know for a fact, that $R\ge r^2$ is always true. Therefore, for me, the answer is:

$$ p = \left\{\begin{array}{1 1}\lfloor (P*r)/R \rfloor & \quad \text{if $P<\lfloor2^{63}/r \rfloor$} \\ \lfloor P/\lfloor R/r\rfloor \rfloor \;[-1] & \quad \text{if $P\ge \lfloor2^{63}/r \rfloor$ and $R\ge r^2$} \\ \text{here be dragons} & \quad \text{if $P\ge \lfloor2^{63}/r \rfloor$ and $R<r^2$} \end{array} \right.\ $$

Here, $[-1]$ indicates that I only subract -1 in some cases, which I check by first leaving it out, testing the value of $p$ I get and then putting it in if I need it.

In the third case, one could use an arbitrary precision approach (c.f. GMP) and calculate the answer by doing long division as suggested by MvG. But maybe someone can find a more elegant way to slay the dragons for completeness... ;) Otherwise, if I have time a little later, I (or someone else) might add the formula for implementing the long division.

Related Solutions

Algorithms – How to Calculate the Number of Decimal Digits for a Binary Number?

I would post this as a comment if I could, because it seems too simple for an answer. You already said that you can find the highest bit with a single instruction (if the number fits into a register). For every result $r$ you get (that is, your number is between $2^r$ and $2^{r+1}-1$), there are only two possible numbers $k$ and $k+1$ of decimal digits. So you could just use $r$ as an index into a lookup table that stores both $k$ and $10^k$, and output $k$ if your number is less than $10^k$, or $k+1$ otherwise.

[Math] How to calculate non-integer exponents

I'll consider the problem of computing $x^\frac1{q}, \; q > 0$; as I've already mentioned in the comments, one can decompose any positive rational number as $m+\dfrac{p}{q}$, where $m,p$ are nonnegative integers, $q$ is a positive integer, and $p < q$. Thus for computing $x^{m+\frac{p}{q}}$, one could use binary exponentiation on $x^m$ and $\left(x^\frac1{q}\right)^p$ and multiply the results accordingly.

A.N. Khovanskiĭ, in his book on continued fractions, displays a continued fraction representation for the binomial function:

$$(1+z)^\alpha=1+\cfrac{2\alpha z}{2+(1-\alpha)z-\cfrac{(1-\alpha^2)z^2}{3(z+2)-\cfrac{(4-\alpha^2)z^2}{5(z+2)-\cfrac{(9-\alpha^2)z^2}{7(z+2)-\cdots}}}}$$

which converges for $|\arg(z+1)| < \pi$.

Letting $z=x-1$ and $\alpha=\dfrac1{q}$, one can then evaluate this continued fraction (with, say, Lentz-Thompson-Barnett) to generate a "seed" that can be subsequently polished with Newton-Raphson, Halley, or any of a number of iterations with high-order convergence. You'll have to experiment with how accurate a seed you need to start up the iteration, by picking a not-too-small tolerance when evaluating the continued fraction.

Here's some Mathematica code demonstrating what I've been saying earlier, for computing $\sqrt[3]{55}$:

With[{q = 3, t = 55, prec = 30},
 y = N[2 + (1 - 1/q) (t - 1), prec];
 c = y; d = 0; k = 1;
 While[True,
  u = (k^2 - q^-2) (t - 1)^2; v = (2 k + 1) (t + 1);
  c = v - u/c; d = 1/(v - u d);
  h = c*d; y *= h;
  If[Abs[h - 1] <= 10^-4, Break[]];
  k++];
 FixedPoint[
  Function[x, x ((1 + q) t - x^q (1 - q))/(x^q (1 + q) - (1 - q) t)], 
  1 + 2 (t - 1)/q/y]]

Here, I've arbitrarily chosen to stop when the continued fraction has already converged to $\approx 4$ digits, and then polished the result with Halley's method. The result here is good to $\approx 28$ digits. Again, you'll have to experiment on the accuracy versus expense of evaluating the "seed", as well as picking the appropriate iteration method for polishing the seed.

Best Answer

Related Solutions

Algorithms – How to Calculate the Number of Decimal Digits for a Binary Number?

[Math] How to calculate non-integer exponents

Related Question