1
$\begingroup$

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

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
10
  • $\begingroup$ Thanks. That is a good suggestion, although right now I do not see any obvious way to solve it. $\endgroup$ – Fabio Nov 4 '16 at 0:39
  • $\begingroup$ I agree, there does not appear to be an obvious solution, but that is what makes it an interesting problem. $\endgroup$ – K. Miller Nov 4 '16 at 2:09
  • $\begingroup$ 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. $\endgroup$ – Fabio Nov 4 '16 at 12:11
  • 1
    $\begingroup$ 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. $\endgroup$ – K. Miller Nov 4 '16 at 15:02
  • $\begingroup$ 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. $\endgroup$ – Fabio Nov 7 '16 at 2:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.