# Solving strict inequalities with integer maths accounting for rounding errors

I have an array of floating point numbers $X_i$ of type T, where T could be either float or double.

These numbers are strictly positive and sorted, i.e.

$$0 < X_0 < X_1 < X_2 < ... < X_n$$

I want to find the smallest floating point number $H$ of type T such that the following inequality is strictly satisfied:

$$\lfloor X_{i+1}H \rfloor > \lfloor X_iH \rfloor \quad i=0\dots n-1$$

Let's assume that the numbers $X_i$ are such that there are no underflow or overflow issues to worry about.

I am approaching the problem as described below, which I think is not robust and, when it works, produces a sub optimal result.

Is there a robust and accurate way? I could accept a sub optimal solution, but at least I need it to be robust.

My current method

I use the notation $m(x)$ to indicate a floating point number which approximates $x$ and may be affected by roundoff error. In the following passages, I modify the inequality below taking as appropriate upper and lower bounds.

Let $a(x)$ be the floating point number closest to the floating point number $x$ toward $-\infty$

Let $b(x)$ be the floating point number closest to the floating point number $x$ toward $+\infty$

I mark with $m()$ the operation affected by rounding error, the original inequality is: $$\lfloor m(X_{i+1}\,H) \rfloor > \lfloor m(X_i\,H )\rfloor ,\quad i=0\dots n-1$$

I remove the $floor$ operator subtracting and adding 0.5. Since this is a conceptual operation, not actually executed on the machine, I do not mark it with an $m()$ $$m(X_{i+1}\,H) - 0.5 > m(X_i\,H) + 0.5 ,\quad i=0\dots n-1$$

I reorganize terms $$m(X_{i+1}\, H) - m(X_i\, H) > + 1 ,\quad i=0\dots n-1$$

I resolve the $m()$ operator with strict LHS minorant and RHS majorant. Note I cannot use $a(x)$ and $b(x)$ here because I do notbexplicitly compute the product $X_i\,H$. I use multiplication by $(1+2\epsilon)$ or $(1-2\epsilon)$ as I believe this will result in a number strictly greater (or lower) than the original one. Does this work? This may be the weak point in my formula! If $X_i$ is very small, multiplying it by $(1+2\epsilon)$ won't change it. Perhaps I should simply take $b(X_i)$ and $a(X_i)$.

$$X_{i+1}\,H\,(1-2\epsilon) - X_{i}\,H\,(1+2\epsilon) > b(1) ,\quad i=0\dots n-1$$

I solve the inequality for $H$, marking the operations affected by rounding error with $m()$

$$H > m\left( \frac{b(1)}{ m( X_{i+1}\,(1-2\epsilon) - X_{i}\,(1+2\epsilon) )} \right) ,\quad i=0\dots n-1$$

Take majorant or minorant as appropriate for rounded terms $$H > b\left( \frac{b(1)} { a( X_{i+1}(1-2\epsilon) - X_{i}(1+2\epsilon) )} \right) ,\quad i=0\dots n-1$$

and eventually take a majorant which gives me the final formula for $H$ $$H = b\left( b\left( \frac{b(1)}{ \underset{0\leq i < n-1}{min}\{ a( X_{i+1}\,(1-2\epsilon) - X_{i}\,(1+2\epsilon) ) \} } \right) \right)$$

I would first try and solve this problem in exact arithmetic before considering floating point arithmetic. This problem is not quite as simple as it seems. For example, consider the sequence $0.5 < 0.7 < 1.1$. According to the exact arithmetic version of your formula, $H = 1/0.2 = 5$. Indeed, this value of $H$ works since $2 = \lfloor 2.5\rfloor < 3 = \lfloor 3.5\rfloor < 5 = \lfloor 5.5 \rfloor$. However, this is not the best value of $H$. For example, $H = 3$ also works, and even smaller values of $H$ are possible. Phrased another way, I believe the question you are trying to answer is what is the smallest value of $H$ such that there exist integers $I_1 < I_2 < I_3 < \ldots < I_n$ with the property that

$$0 < HX_0 < I_1 \leq HX_1 < I_2 \leq HX_2 < I_3 \leq \ldots \leq HX_{n-1} < I_n \leq HX_n$$

Since $I_n \geq n$, we get that $H$ is bounded below by $n/X_n$. Thus, for any $H < n/X_n$, you are guaranteed that $\lfloor HX_{i+1} \rfloor = \lfloor HX_i\rfloor$ for some $i = 0,\ldots,n-1$. Therefore, you have the following bounds

$$\frac{n}{X_n} \leq H \leq \max_{i=0,\ldots,n-1}\left(\frac{1}{X_{i+1} - X_i}\right)$$

As Fabio pointed out this lower bound can be tightened to

$$\max_{i=1,\ldots,n}\left(\frac{i}{X_i}\right) \leq H$$

With this lower bound, I propose the following simple algorithm to compute the optimal value of $H$.

Let x(1),...,x(n) be positive numbers ordered
so that x(1) < x(2) < ... < x(n)
L = max(2/x(2),...,n/x(n))
H = L
tol = 1 + eps
nsteps = 0
while (true)
nsteps = nsteps + 1
k = 0
for i = 2,...,n do
if floor(H * x(i)) <= floor(H * x(i - 1)) then
k = i
break
end
end
if k equals 0 then
break
end
H = ceiling(H * x(k)) / x(k)
while floor(H * x(k)) equals floor(H * x(k - 1)) do
H = H * tol
end
end


Multiplying $H$ by the tolerance $\text{tol} = 1 + \epsilon$, where $\epsilon$ is machine epsilon, accounts for any roundoff error that makes $H$ slightly too small. The cost of this algorithm is $O(nM)$, where $M$ is the number of steps. Since $H$ is strictly increasing the algorithm must terminate in a finite number of steps, i.e., when $H$ reaches the upper bound. However, finding a reasonable bound for $M$ seems difficult as it appears to be highly dependent on the distribution of points $\{X_i\}$ as well as the number of points $n$. Some preliminary testing suggests that this algorithm finds the best $H$, however, at this point I do not have a proof that it gives the optimal value.

• Thanks. That is a good suggestion, although right now I do not see any obvious way to solve it. – Fabio Nov 4 '16 at 0:39
• I agree, there does not appear to be an obvious solution, but that is what makes it an interesting problem. – K. Miller Nov 4 '16 at 2:09
• I think the lower bound could be tightened to max( i/X[i] ), because it must be I[n]>n for any i. However this can originate pairs HX[i], HX[i+1] which are less than 1 unit apart and resolve to different integers purely because they are slightly on the left and of the right of an integer. Here a tiny rounding error can make them move fall on the same side, causing problems. – Fabio Nov 4 '16 at 12:11
• Good call, the lower bound can definitely be improved. I will update my answer to reflect this. I have also devised a simple algorithm that I believe will calculate the minimal feasible $H$. I will post it shortly. – K. Miller Nov 4 '16 at 15:02
• Thanks, I had the same thought. Only suggestion is to replace $H=H\,tol$ with $H=next(H)$. In fact, for small $H$, it can be $H\,tol \equiv tol$. Take for instance $H=m$, where $m$ is the smallest positive floating point number, then $H(1+eps) \equiv H$. For my use case however, I would prefer something more predictable in terms of running time. I came up with a robust and simpler version of the algorithm I posted previously. Still suboptimal, but robust and with complexity close to $O(N)$. I will update my question. – Fabio Nov 7 '16 at 2:50