x(x-1)/ 2 = c的快整数解

Question

给定一个非负整数c,我需要一个高效的算法来找到最大的整数x使得

x*(x-1)/2 <= c

同样,我需要一个高效且可靠的精确算法来计算:

x = floor((1 + sqrt(1 + 8*c))/2)        (1)

为了定义,我标记了这个问题C++,所以答案应该是用该语言编写的函数.您可以假设这c是一个无符号的32位int.

此外,如果您可以证明(1)(或涉及浮点运算的等效表达式)总是给出正确的结果,那也是一个有效的答案,因为现代处理器上的浮点数可能比整数算法更快.

Answer 1

如果你愿意假设IEEE为所有操作(包括平方根)进行正确的舍入,那么你编写的表达式(加上强制转换为double)会在所有输入上给出正确的答案.

这是一个非正式的证明.由于c32位无符号整数被转换为具有53位有效数的浮点类型,因此1 + 8*(double)c是精确的,并且sqrt(1 + 8*(double)c)是正确舍入的.1 + sqrt(1 + 8*(double)c)精确到一个ulp内,因为最后一个项小于2**((32 + 3)/2) = 2**17.5意味着后一项最后一位的单位小于1,因此(1 + sqrt(1 + 8*(double)c))/2精确到一个ulp内,因为除法2是精确的.

最后一项业务是发言.这里的问题(1 + sqrt(1 + 8*(double)c))/2是将舍入为整数的时间.当且仅当sqrt(...)舍入到奇数时才会发生这种情况.由于参数sqrt是一个整数,最坏的情况看起来像sqrt(z**2 - 1)正奇数整数z,我们绑定

z - sqrt(z**2 - 1) = z * (1 - sqrt(1 - 1/z**2)) >= 1/(2*z)

由泰勒扩张.由于z小于2**17.5,与最接近整数的间隙至少1/2**18.5在幅度小于的结果上2**17.5,这意味着该错误不能由正确舍入产生sqrt.

采用Yakk的简化,我们可以写

(uint32_t)(0.5 + sqrt(0.25 + 2.0*c))

没有进一步检查.