Find smallest integer that satisfies floating point inequality equation

Question

I am looking for a fast algorithm that finds the smallest integer N that will satisfy the following inequality where s, q, u, and p are float numbers (using the IEEE-754 binary32 format):

s > q + u * p / (N - 1)

where N can be any positive integer represented by a signed 32-bit integer. After (N - 1) is converted to float, all arithmetic is evaluated in float.

Additional constraints are:

• 0 < p < 1.
• -1 ≤ q ≤ 1.
• q < s.
• 0 < u.

I am having trouble figuring out how to do this in a robust way that deals with floating point rounding errors and comparison properly. Here is my poor attempt at a solution that is not fast and is not even robust in that I cannot determine the minimum SOME_AMOUNT:

int n = std::max(1.0f, floorf((u * p / (s - q)) - 1.0f));

// Floating point math might require to round up by some amount...
for (int i = 0; i < SOME_AMOUNT; ++i)
    if (!(q + (u * p / (n + 1)) < second))
        ++n;

You can see above my formula for calculating n using basic algebra. The for loop is my crude means of trying to account for floating point rounding errors. I am checking it with brute force like this:

int nExact = 0;
bool found = false;
for (; nExact < SOME_BIG_NUMBER; ++nExact) {
    if (q + (u * p / (nExact + 1)) < second) {
        found = true;
        break;
    }
}
assert(found);
assert(n == nExact);

Any floating point gurus have a reasonably fast an answer in C++?

Frankly, if someone can even just give a theoretically sound proof of an upper bound on 'SOME_AMOUNT' above I'd be reasonably happy...

Answer 1

To be on safe side, we can first get a larger possible value (upper bound) and a smaller possible value (lower bound) and than reduce it to our actual answer, this way it will be accurate and faster that just iterating over numbers.

By solving the inequality we get,

N > u * p / (s - q) + 1

Getting an upper bound

So you will first find a maximum guessed answer, by using integers. We will increase numerator and integer cast denominator

int UP = (int)(u * p + 1);    // Increase by one
int D = (int)(s - q);         // we don't increase this because it  would cause g to decrease, which we don't want

float g = UP / (float)D + 1;  // we again float cast D to avoid integer division
int R = (int)(g + 1);         // Now again increase g

/******** Or a more straight forward approach ********/
int R = (int)(((int)(u*p+1))/(s-q) + 1 + 1)

// Add rounding-off error here
if(R + 128 < 0) R = 2147483647;    // The case of overflow
else R += 128;

This is your maximum answer (upper bound).

Getting a lower bound

Just as previous but this time we will increase denominator and integer cast numerator

int UP = (int)(u * p);         // will automatically decrease
int D = (int)(s - q + 1);      // we increase this because it would cause g to decrease, which we want

float g = UP / (float)D + 1;   // we again float cast D to avoid integer division
int L = (int)g;                // Integer cast, will automatically decrease
/******** Or a more straight forward approach ********/
int L = (int)(((int)(u*p))/(s-q+1) + 1)

// Subtract rounding-off error
if(L - 128 <= 1 ) L = 2;        // N cannot be below 2
else L -= 128;

This is your minimum answer (lower bound).

Note: The reason of integer casting is to reduce our sample space. It can be omitted if you feel so.

Elimination of possible numbers and getting the right one

for (int i = L; i <= R; ++i){
    if ((s > q + u*p/(i-1))) break;   // answer would be i
}
N = i;    // least number which satisfies the condition

You can do it even faster with binary search if the gap between bounds (R-L) is large. As for number-range whose difference is 2^n can be reduced in only n steps.

// we know that
// lower limit = L;
// upper limit = R;
// Declare u, p, q, s in global space or pass as parameters to biranySearch

int binarySearch(int l, int r)
{
    if(l==r) return l;

    if (r > l) {
        int mid = l + (r - l) / 2;

        bool b = (s > q + (p*u)/(mid-1));

        if (b==true){
            // we know that numbers >= mid will all satisfy
            // so our scope reduced to [l, mid]
            return binarySearch(l, mid);
        }
        // If mid doesn't satisfy
        // we know that our element is greater than mid
        return binarySearch(mid+1, r); 
    } 
} 

int main(void) 
{
    // calculate lower bound L and upper bound R here using above methods
    int N = binarySearch(L, R);
    // N might have rounding-off errors, so check for them
    // There might be fluctuation of 128 [-63 to 64] so we will manually check.
    // To be on safe side I will assume fluctuation of 256
    L = N-128 > 2 ? N-128 : 2;
    R = N+128 < 0 ? 2147483647 : N+128;
    for(int i=L; i<=R; ++i){
        if( s > q + u * p / ((float)i - 1)) {
            break;
        }
    }
    cout << i << endl;
}

It is mostly a concept, but it is fast and safe. The only thing is that I have not tested it, but it should work!

Answer 2

Here is the start of a solution. Some caveats:

• It is in C, not C++.
• It assumes IEEE-754 arithmetic with rounding to nearest.
• It does not handle cases where the inequality requires N to go out of the bounds from 2 to INT_MAX.
• I have not tested it much.

The code first using floating-point arithmetic to estimate where the boundary where the inequality changes is, neglecting rounding errors. It tests the inequality to see whether it needs to increase or decrease the candidate value. Then it iterates through consecutive integer float values to find the boundary. My sense is this will take few iterations, but I have not analyzed it completely.

This produces the least float with an integer value that satisfies the inequality when used in place of the denominator N-1. The code then finds the least int N such that N-1 rounds to that float, and that should be the N that is the least int for which the inequality is satisfied.

#include <math.h>
#include <stdio.h>
#include <stdlib.h>


//  Test the inequality.
static int Test(float s, float q, float u, float p, int N)
{
    return s > q + (float) (((float) (u * p)) / (N-1));
}


int main(void)
{
    float s = 1;
    float q = 0;
    float u = 0x1p30, p = 1;

    /*  Approximate the desired denominator (N-1) -- would be exact with real
        arithmetic but is subject to rounding errors.
    */
    float D = floorf(u*p/(s-q));

    //  Test which side of the boundary where the inequality changes we are on.
    if (Test(s, q, u, p, (int) D + 1))
    {
        //  We are above the boundary, decrement find the boundary.
        float NextD = D;
        do
        {
            D = NextD;
            //  Decrement D by the greater of 1 or 1 ULP.
            NextD = fminf(D-1, nexttowardf(D, 0));
        }
        while (Test(s, q, u, p, (int) NextD + 1));
    }
    else
        //  We are below the boundary, increment to find the boundary.
        do
            //  Increment D by the greater of 1 or 1 ULP.
            D = fmaxf(D+1, nexttowardf(D, INFINITY));
        while (!Test(s, q, u, p, (int) D + 1));

    //  Find the distance to the next lower float, as an integer.
    int distance = D - nexttowardf(D, 0);

    /*  Find the least integer that rounds to D.  If the distance to the next
        lower float is less than 1, then D is that integer.  Otherwise, we want
        either the midpoint between the D and the next lower float or one more
        than that, depending on whether the low bit of D in the float
        significand is even (midpoint will round to it, so use midpoint) or odd
        (midpoint will not round to it, so use one higher).

        (int) D - distance/2 is the midpoint.

        ((int) D / distance) & 1 scales D to bring the low bit of its
        significand to the one’s position and tests it, producing 0 if it is
        even and 1 if it is odd.
    */
    int I = distance == 0 ? (int) D
        : (int) D - distance/2 + (((int) D / distance) & 1);

    //  Set N to one more than that integer.
    int N = I+1;

    printf("N = %d.\n", N);

    if (Test(s, q, u, p, N-1) || !Test(s, q, u, p, N))
    {
        fprintf(stderr, "Error, solution is wrong.\n");
        exit(EXIT_FAILURE);
    }
}