给定 (a, b) 计算 k 的最大值，使得 a^{1/k} 和 b^{1/k} 是整数答案

【问题标题】：Given (a, b) compute the maximum value of k such that a^{1/k} and b^{1/k} are whole numbers给定 (a, b) 计算 k 的最大值，使得 a^{1/k} 和 b^{1/k} 是整数
【发布时间】：2019-01-02 15:16:46
【问题描述】：

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

这是我的代码的 sn-p，我已对其进行了注释以进行澄清。

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 个根都等于 0 mod 1（这意味着它们可以被 1 整除，因此整数）。

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

我尝试过的另一种方法是更改循环条件以检查 a^{1/k} 的下限是否等于 a^{1/k} 本身。但同样，这会无限运行，可能是由于精度错误。

有人知道我该如何解决这个问题吗？

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

【问题讨论】：

对于浮点数，测试相等性是不合理的。见floating-point-gui.de
为什么一个返回值针对1.0 进行测试，而另一个只针对0？
阅读ericlippert.com/2014/03/05/how-to-debug-small-programs 了解调试代码的技巧。
k的最小值，还是1/k的最小值（对应k的最大值）？后者似乎更有趣，因此更有可能成为要解决的实际问题。
（我看到你删除了你的question about the two integer sum。为了它的价值，我试图修复我的代码来处理你找到的反例。这里：ideone.com/JKRkqy）

标签： c++ algorithm c++11 math

【解决方案1】：

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

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

所以最快的方法可能是首先检查a 和b 是否为整数，然后执行prime decomposition 使得a=p₀^e₀ * p₁^e₁ * ...，其中 p_i 是不同的素数。

请注意，要使 a^1/k 成为整数，每个 e_i 也必须能被 k 整除。换句话说，k 必须是 e_i 的公约数。如果 b^1/k 是整数，则 b 的素数幂也必须如此。

因此，最大的k 是a 和b 的所有e_i 中的greatest common divisor。

使用您的方法，您将遇到大量问题。所有 IIEEE 754 binary64 浮点（x86 上双精度的情况）都有 53 个有效位。这意味着所有大于 2⁵³ 的 double 都是整数。

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

另一方面，由于只有 53 个有效位，所以 double 的素数分解是微不足道的。

【讨论】：

【解决方案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) 函数（并使用它而不是 ==，并使用它的否定而不是 !=）。

根据经验，永远不要测试浮点数是否相等（即==），而是考虑它们之间足够小的距离；也就是说，将x == y（分别为x != y）之类的测试替换为fabs(x-y) < EPSILON（分别为fabs(x-y) > EPSILON）之类的东西，其中EPSILON是一个小的正数，因此测试一个小的L₁ distance（为了相等，和足够大的距离以防止不平等）。

并避免整数问题中的浮点数。

实际上，预测或估计浮点精度非常困难。您可能需要考虑像CADNA 这样的工具。我的同事 Franck Védrine 是估计数值误差的静态程序分析器方面的专家（参见例如他的 TERATEC 2017 presentation on Fluctuat）。这是一个困难的研究课题，另见 D.Monniaux 的论文the pitfalls of verifying floating-point computations 等。

浮点错误did in some cases 会造成人命损失（或损失数十亿美元）。搜索网络以获取详细信息。在某些情况下，计算出的数字所有数字都是错误的（因为错误可能会累积，最终的结果是通过数千次运算得到的）！和chaos theory有一些间接关系，因为很多程序可能有一些numerical instability。

【讨论】：

好的，知道了。但只有一件事——循环条件不应该是“fabs(fmod(pow(a, 1.0/k), 1) - 1.0) >= 0.0000001”，因为我们想在语句不正确时继续循环吗？
也是绝对值，我们应该减去 0.0 而不是 1.0 对吧？因为它们与 0 mod 1 一致？只是确保...
这是关于计算距离的。您甚至可以使用euclidean distance，但计算时间会更长。
是的。目前，布尔值测试 fmod(pow(a, 1.0/k), 1) 的值是否在 1 的 epsilon 范围内。但它不应该测试它是否在 0 的 epsilon 范围内吗？因为我们希望 while 条件在它们都是整数时评估为“假”，这恰好发生在 a 和 b 的第 k 个根等于 0 mod 1，而不是 1 mod 1
您希望fmod(pow(a, 1.0/k), 1)（称为x）与1.0（称为y）有足够的不同。所以你希望fabs(x-y) 大于一些EPSILON

【解决方案3】：

正如其他人所提到的，比较浮点值是否相等是有问题的。如果你找到一种直接处理整数的方法，你就可以避免这个问题。一种方法是将整数提高到k 的幂，而不是取kth 的根。细节留给读者练习。

【讨论】：