给定（a，b），计算最大的k值，使得a^{1/k}和b^{1/k}都是整数。

Question

我正在编写一个程序，尝试找到最小的k值（k > 1），使得a和b（两者都已给定）的k次方根等于整数。

这是我的代码片段，我已经添加了注释以便更好地理解：

int main()
{
    // Declare the variables a and b.
    double a;
    double b;
    // Read in variables a and b.
    while (cin >> a >> b) {

        int k = 2;

        // We require the kth root of a and b to both be whole numbers.
        // "while a^{1/k} and b^{1/k} are not both whole numbers..."
        while ((fmod(pow(a, 1.0/k), 1) != 1.0) || (fmod(pow(b, 1.0/k), 1) != 0)) {

        k++;

        }

基本上，我读入（a，b），并从k = 2开始递增k，直到a和b的k次根都模1同余（这意味着它们可以被1整除并且是整数）。

但是，循环运行无限。我尝试过研究，我认为这可能与精度误差有关，但我不太确定。

我尝试的另一种方法是将循环条件更改为检查a ^ {1 / k}的底部是否等于a ^ {1 / k}本身。但是，这再次无限运行，可能是由于精度误差所致。

有人知道我怎样才能解决这个问题吗？

编辑：例如，当（a，b）=（216，125）时，我希望k = 3，因为216 ^（1/3）和125 ^（1/3）都是整数（即5和6）。

- oompa

---

使用浮点数进行等值测试是不合理的。参见http://floating-point-gui.de/。 - Basile Starynkevitch

为什么一个返回值要与 1.0 进行比较，而另一个只与 0 进行比较？ - Tim Randall

阅读 https://ericlippert.com/2014/03/05/how-to-debug-small-programs/ 了解有关调试代码的技巧。 - Code-Apprentice

k的最小值，还是1/k的最小值（对应于k的最大值）？后者似乎更有趣，因此更可能是要解决的实际问题。 - President James K. Polk

我看到你删除了关于两个整数求和的问题（https://stackoverflow.com/questions/58150851/finding-two-integers-without-a-zero-digit-that-sum-to-a-target-integer）。值得一提的是，我尝试修复我的代码以处理你发现的反例。这里是链接：https://ideone.com/JKRkqy - גלעד ברקן

Answer 1

这不是一个编程问题，而是一个数学问题：

如果a是实数，k是正整数，并且a^(1./k)是整数，则a是整数。（否则目的是为了玩弄近似误差）

因此，最快的方法可能是先检查a和b是否为整数，然后进行质因数分解，使a=p0^e0 * p1^e1 * ...，其中p_i是不同的质数。

注意，为了a^（1/k）是整数，每个ei也必须被k整除。换句话说，k必须是ei的公共因子。如果b^（1/k）也是整数，则对b的质数幂也必须如此。

因此，最大的k是a和b的所有e_i的最大公约数。

---

如果您采用这种方法，处理大数字时可能会遇到问题。所有IIEEE 754二进制64位浮点数(double在x86上的情况)都有53个有效位。这意味着大于2⁵³的所有双精度浮点数都是整数。

pow(x,1./k)函数对于两个不同的x产生相同的值，因此使用您的方法将会得出错误的答案。例如，数字5⁵*2⁹⁰和3⁵*2¹²⁰可以用double精确表示。该算法的结果为k=5。您可能会在这些数字中找到k=5，但是您也会在5⁵*2⁹⁰-2⁴⁹和3⁵*2¹²⁰中找到k=5，因为pow(5⁵*2⁹⁰-2⁴⁹,1./5)==pow(5⁵*2⁹⁰)。演示在这里。

另一方面，由于只有53位有效位数，因此双精度数的质因数分解是平凡的。

- Oliv

Answer 2

浮点数不是数学上的实数。计算是“近似”的。请参见http://floating-point-gui.de/ 您可以将测试fmod(pow(a, 1.0/k), 1) != 1.0替换为类似于fabs(fmod(pow(a, 1.0/k), 1) - 1.0) > 0.0000001的内容（并尝试使用各种这样的值，而不是0.0000001；另请参见std::numeric_limits::epsilon，但要谨慎使用，因为pow在计算中可能会产生一些误差，而1.0/k也会注入不精确性——详细信息非常复杂，请深入IEEE754规范）。

当然，您可以（并且可能应该）定义您的bool almost_equal(double x, double y)函数（并使用它代替==，并使用其否定代替!=）。

作为一个经验法则，永远不要测试浮点数是否相等（即==），而是考虑它们之间足够小的距离；也就是说，用类似于fabs(x-y) < EPSILON（或者fabs(x-y) > EPSILON）来替换像x == y（或者x != y）这样的测试，其中EPSILON是一个小正数，因此测试一个小的L₁ distance（用于相等性），以及一个足够大的距离用于不等式。

并且避免在整数问题中使用浮点数。

实际上，预测或估计浮点精度非常困难。您可能需要考虑使用像CADNA这样的工具。我的同事Franck Védrine是静态程序分析器估计数值误差方面的专家（例如，请参见他在TERATEC 2017 Fluctuat演示文稿中的介绍）。这是一个困难的研究课题，也请参阅D.Monniaux的论文验证浮点计算的陷阱等。

浮点误差有时会导致人员伤亡（或数十亿美元的损失）。请在网上搜索详细信息。有些情况下，计算得到的数字中所有位数都是错误的（因为误差可能会累积，并且最终结果是通过组合成千上万个操作得到的！）这与混沌理论有一些间接关系，因为许多程序可能存在数值不稳定性。

- Basile Starynkevitch

---

好的，明白了。但是只有一件事——循环条件应该是“fabs(fmod(pow(a, 1.0/k), 1) - 1.0) >= 0.0000001”，因为我们希望在语句不成立时继续循环。 - oompa

在绝对值中，我们应该减去0.0而不是1.0，对吧？因为它们对于mod 1同余于0？只是确认一下... - oompa

这是关于计算距离的。你甚至可以使用欧几里得距离，但这需要更长的计算时间。 - Basile Starynkevitch

是的。目前，布尔值测试fmod(pow(a, 1.0/k), 1)的值是否在epsilon范围内等于1。但是，它不应该测试它是否在epsilon范围内等于0吗？因为我们希望当a和b的k次方根等于0 mod 1而不是1 mod 1时，while条件求值为“false”，这种情况恰好发生在它们都是整数时。 - oompa

你想让 fmod(pow(a, 1.0/k), 1)（称之为 x）与 1.0（称之为 y）有足够的差异。因此，你希望 fabs(x-y) 大于某个 EPSILON。 - Basile Starynkevitch

Answer 3

正如其他人所提到的，比较浮点数的相等性是有问题的。如果你找到一种直接使用整数的方法，就可以避免这个问题。其中一种方法是将整数提高到k次幂，而不是取k次根。具体细节留给读者自行练习。

- Code-Apprentice