Python中的费马大定理算法

Question

我在 Python 中创建了一个程序来寻找 Fermat 大定理的解决方案（我知道这无法解决，但我只是将它用作编程练习）。费马大定理指出：

对于任何大于 2 的 n 整数值，没有三个正整数 a、b 和 c 满足方程 a^n + b^n = c^n。

来源：Wikipedia

我的算法（在 Python 中）在这里：

from fractions import Fraction

def root(num):
    return num ** Fraction(1 / power)

def two_numbers():
    a = b = 1
    while True:
        yield a, b
        if a == b:
            a += 1
            b = 1
        else:
            b += 1

def test(a, b):
    return root((a ** power) + (b ** power)) % 1 == 0

power = int(input('Power: '))

generator = two_numbers()
for a, b in generator:
    result = test(a, b)
    print(a, b, result)
    if result:
        break

当使用 12 的幂（在提示符处输入 12）运行时，它会停在：

17 1 True

这可能是由于小数索引不准确。

我怎么知道它是否被四舍五入，或者以其他方式解决这个问题？

Answer 1

Fraction 不能表示无理数，所以即使使用也会遇到舍入错误。完全避免浮点数不会给舍入错误留下空间：

def is_perfect_kth_power(n, k):
    low = 1
    high = 1

    # Find an upper bound for the binary search
    while high**k < n:
        high *= 2

    while low + 1 < high:
        midpoint = (low + high) // 2

        if n < midpoint**k:
            high = midpoint
        elif midpoint**k < n:
            low = midpoint
        else:
            return True

    return False

由于您试图找到不存在的东西，我认为运行时不是这里的问题。

Answer 2

是的，相对误差的一阶项为 17

(17**12+1)**(1.0/12)

是

1.0/(p*a^p) = 1.0/(12*17**12) = 1.43031501388558e-16. 

这小于机器 epsilon 2.2e-16，即，不足以影响分数幂计算中使用的 64 位浮点数的尾数。

---

你可以使用

def test(a, b):
    num=a**power+b**power; 
    c=root(num); 
    return num - int(c+0.5)**power == 0

直到由于从浮点数到整数的转换中的溢出问题而失败。

---

您可以反转一阶相对误差的计算。你从等式中得到c 大于a，因此作为整数你需要c 至少a+1。插入给出

 a**p+b**p >= (a+1)**p = a**p + p*a**(p-1) + ...

由二项式定理，这样你至少想要

b**p > p*a**(p-1)  <==>  b > a*(p/a)**(1/p)

使用这个下限应该可以避免这些舍入问题。

对于a=17，这将17 作为b 的下限，因此在b<a 规则下没有尝试的情况。对于a=171，下限是138，实际上c>=172 的第一种情况是b=138。

Answer 3

我会使用以下策略来继续使用整数：

• 对于所需范围内的所有a 值，计算a^n，

• 扫描b在所需范围内的所有值并计算a^n + b^n；同时，维护一个变量c，并确保始终保持(c-1)^n < a^n + b^n <= c^n。

随着a^n + b^n 的增加，只有右手不等式会失效，您可以根据需要多次增加c 来解决此问题。要初始化c，只需注意a^n < a^n + 1 并以c= a 开头即可。

n= 2
m= 50
for a in range(1, m+1):
    an= a ** n
    c= a; cn= an
    for b in range(a+1, m+1):
        anbn= an + b ** n
        while anbn > cn:
            c+= 1
            cn= c ** n
        if anbn == cn: # Bingo!
            print a, b, c

对于n=2的情况，

3 4 5
5 12 13
6 8 10
7 24 25
8 15 17
9 12 15
9 40 41
10 24 26
12 16 20
12 35 37
14 48 50
15 20 25
15 36 39
16 30 34
18 24 30
20 21 29
20 48 52
21 28 35
24 32 40
24 45 51
27 36 45
28 45 53
30 40 50
33 44 55
36 48 60
40 42 58
...

（请注意，我们从b= a+1 开始以避免重复。）

Answer 4

在这里，您的问题来自精确度。您可以通过以下方式检查：

print(np.finfo(float).eps)
# 2.22044604925e-16

或

print(np.finfo(np.float32).eps)
# 1.19209e-07

为了更精确，您可以使用“十进制”库并更新您的根函数：

from decimal import Decimal

def root(num, power):
    return Decimal(num ** Decimal(1 / power))