我正在寻找一种快速算法，该算法可以找到满足以下不等式的最小整数N，其中s，q，u和p是float数字(使用IEEE-754二进制32格式(：
s > q + u * p / (N - 1)
其中 N 可以是由有符号 32 位整数表示的任何正整数。(N - 1)转换为float后，所有的算术都用float计算。
其他约束包括：
- 0<
p><1.> - -1 ≤
q≤ 1. q- 0
我无法弄清楚如何以稳健的方式正确处理浮点舍入误差和比较。 这是我对解决方案的糟糕尝试，该解决方案速度不快，甚至不健壮，因为我无法确定最小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"的上限，我会相当高兴......
这是一个解决方案的开始。一些注意事项：
- 它是 C 语言，而不是 C++。
- 它假定 IEEE-754 算术，四舍五入到最接近。
- 它不处理不等式要求 N 超出从 2 到
INT_MAX的范围的情况。 - 我没有测试太多。
代码首先使用浮点运算来估计不等式变化的边界在哪里，忽略舍入误差。它测试不等式以查看是否需要增加或减少候选值。然后，它循环访问连续的整数float值以查找边界。我的感觉是这将需要几次迭代，但我还没有完全分析它。
这将产生最小float，其整数值在代替分母N-1时满足不等式。然后代码找到最小int N，使得N-1舍入到该float，这应该是满足不等式的最不intN。
#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);
}
}
为了安全起见，我们可以首先获得一个更大的可能值(上限(和一个更小的可能值(下限(，然后将其减少到我们的实际答案，这样它将比仅迭代数字更准确、更快。
通过解决我们得到的不平等，
N > u * p / (s - q) + 1
获取上限
因此，您将首先通过使用整数找到最大猜测答案。我们将增加分子和整数铸造分母
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;
这是您的最大答案(上限(。
获得下限
和以前一样，但这次我们将增加分母和整数铸造分子
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;
这是您的最低答案(下限(。
注意：整数转换的原因是减少我们的样本空间。如果你觉得这样，可以省略。
消除可能的数字并获得正确的数字
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
如果边界之间的间隙 (R-L( 很大，您可以使用二叉搜索更快地完成此操作。至于相差为 2^n 的数字范围只能用 n 步减少。
// 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;
}
它主要是一个概念，但它既快速又安全。唯一的问题是我没有测试过它，但它应该可以工作！