c++ - 找到满足浮点不等式的最小整数

Question

我正在寻找一种快速算法来找到满足以下不等式的最小整数 N，其中 s、q、u、和 p 是 float 数字(使用 IEEE-754 binary32 格式):

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

其中 N 可以是由带符号的 32 位整数表示的任何正整数。 (N - 1) 转换为 float 后，所有算术运算都在 float 中计算。

附加约束是:

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

我无法弄清楚如何以可靠的方式执行此操作，以正确处理浮点舍入错误和比较。这是我对解决方案的拙劣尝试，它速度不快，甚至不够稳健，因为我无法确定最小 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;

您可以在上面看到我使用基本代数计算 n 的公式。 for 循环是我尝试解决浮点舍入错误的粗略方法。我正在用这样的蛮力检查它:

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

有哪位浮点大师在 C++ 中有相当快的答案？

坦率地说，如果有人甚至可以就上面的“SOME_AMOUNT”上限提供理论上可靠的证据，我会相当高兴......

Answer 1

这是一个解决方案的开始。一些注意事项:

• 它是用 C 语言编写的，而不是 C++。
• 它采用四舍五入到最接近的 IEEE-754 算法。
• 它不处理不等式要求 N 超出范围从 2 到 INT_MAX 的情况。
• 我没有对其进行太多测试。

代码首先使用浮点运算来估计不等式变化的边界在哪里，忽略舍入误差。它测试不等式以查看是否需要增加或减少候选值。然后它遍历连续的整数 float 值以找到边界。我的感觉是这将需要几次迭代，但我还没有完全分析它。

当用于代替分母 N-1 时，这会产生最小的 float，其整数值满足不等式。然后代码找到最小的 int N 使得 N-1 舍入到那个 float，那应该是 N 是满足不等式的最小 int。

#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);
    }
}