Solving strict inequalities with integer maths accounting for rounding errors

Question

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)$$

Answer 1

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.