forked from norvig/pytudes
-
Notifications
You must be signed in to change notification settings - Fork 0
/
beal.py
159 lines (135 loc) · 6.03 KB
/
beal.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
"""Search for counterexamples to Beal's conjecture
See http://norvig.com/beal.html and http://www.bealconjecture.com"""
from __future__ import division, print_function
from math import log
from itertools import combinations, product
from collections import defaultdict
try:
from math import gcd # For Python 3.6 and up
except ImportError:
from fractions import gcd # For older versions (works in 2.7 as well)
def beal(max_A, max_x):
"""See if any A ** x + B ** y equals some C ** z, with gcd(A, B) == 1.
Consider any 1 <= A,B <= max_A and x,y <= max_x, with x,y prime or 4."""
Apowers = make_Apowers(max_A, max_x)
Czroots = make_Czroots(Apowers)
for (A, B) in combinations(Apowers, 2):
if gcd(A, B) == 1:
for (Ax, By) in product(Apowers[A], Apowers[B]):
Cz = Ax + By
if Cz in Czroots:
C = Czroots[Cz]
x, y, z = exponent(Ax, A), exponent(By, B), exponent(Cz, C)
print('{} ** {} + {} ** {} == {} ** {} == {}'
.format(A, x, B, y, C, z, C ** z))
def make_Apowers(max_A, max_x):
"A dict of {A: [A**3, A**4, ...], ...}."
exponents = exponents_upto(max_x)
return {A: [A ** x for x in (exponents if (A != 1) else [3])]
for A in range(1, max_A+1)}
def make_Czroots(Apowers): return {Cz: C for C in Apowers for Cz in Apowers[C]}
def exponents_upto(max_x):
"Return all odd primes up to max_x, as well as 4."
exponents = [3, 4] if max_x >= 4 else [3] if max_x == 3 else []
for x in range(5, max_x, 2):
if not any(x % p == 0 for p in exponents):
exponents.append(x)
return exponents
def exponent(Cz, C):
"""Recover z such that C ** z == Cz (or equivalently z = log Cz base C).
For exponent(1, 1), arbitrarily choose to return 3."""
return 3 if (Cz == C == 1) else int(round(log(Cz, C)))
##############################################################################
def tests():
assert make_Apowers(6, 10) == {
1: [1],
2: [8, 16, 32, 128],
3: [27, 81, 243, 2187],
4: [64, 256, 1024, 16384],
5: [125, 625, 3125, 78125],
6: [216, 1296, 7776, 279936]}
assert make_Czroots(make_Apowers(5, 8)) == {
1: 1, 8: 2, 16: 2, 27: 3, 32: 2, 64: 4, 81: 3,
125: 5, 128: 2, 243: 3, 256: 4, 625: 5, 1024: 4,
2187: 3, 3125: 5, 16384: 4, 78125: 5}
Czroots = make_Czroots(make_Apowers(100, 100))
assert 3 ** 3 + 6 ** 3 in Czroots
assert 99 ** 97 in Czroots
assert 101 ** 100 not in Czroots
assert Czroots[99 ** 97] == 99
assert exponent(10 ** 5, 10) == 5
assert exponent(7 ** 3, 7) == 3
assert exponent(1234 ** 999, 1234) == 999
assert exponent(12345 ** 6789, 12345) == 6789
assert exponent(3 ** 10000, 3) == 10000
assert exponent(1, 1) == 3
assert exponents_upto(2) == []
assert exponents_upto(3) == [3]
assert exponents_upto(4) == [3, 4]
assert exponents_upto(40) == [3, 4, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
assert exponents_upto(100) == [
3, 4, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61,
67, 71, 73, 79, 83, 89, 97]
assert gcd(3, 6) == 3
assert gcd(3, 7) == 1
assert gcd(861591083269373931, 94815872265407) == 97
assert gcd(2*3*5*(7**10)*(11**12), 3*(7**5)*(11**13)*17) == 3*(7**5)*(11**12)
return 'tests pass'
##############################################################################
def beal_modp(max_A, max_x, p=2**31-1):
"""See if any A ** x + B ** y equals some C ** z (mod p), with gcd(A, B) == 1.
If so, verify that the equation works without the (mod p).
Consider any 1 <= A,B <= max_A and x,y <= max_x, with x,y prime or 4."""
assert p >= max_A
Apowers = make_Apowers_modp(max_A, max_x, p)
Czroots = make_Czroots_modp(Apowers)
for (A, B) in combinations(Apowers, 2):
if gcd(A, B) == 1:
for (Axp, x), (Byp, y) in product(Apowers[A], Apowers[B]):
Czp = (Axp + Byp) % p
if Czp in Czroots:
lhs = A ** x + B ** y
for (C, z) in Czroots[Czp]:
if lhs == C ** z:
print('{} ** {} + {} ** {} == {} ** {} == {}'
.format(A, x, B, y, C, z, C ** z))
def make_Apowers_modp(max_A, max_x, p):
"A dict of {A: [(A**3 (mod p), 3), (A**4 (mod p), 4), ...]}."
exponents = exponents_upto(max_x)
return {A: [(pow(A, x, p), x) for x in (exponents if (A != 1) else [3])]
for A in range(1, max_A+1)}
def make_Czroots_modp(Apowers):
"A dict of {C**z (mod p): [(C, z),...]}"
Czroots = defaultdict(list)
for A in Apowers:
for (Axp, x) in Apowers[A]:
Czroots[Axp].append((A, x))
return Czroots
##############################################################################
def simpsons(bases, powers):
"""Find the integers (A, B, C, n) that come closest to solving
Fermat's equation, A ** n + B ** n == C ** n.
Let A, B range over all pairs of bases and n over all powers."""
equations = ((A, B, iroot(A ** n + B ** n, n), n)
for A, B in combinations(bases, 2)
for n in powers)
return min(equations, key=relative_error)
def iroot(i, n):
"The integer closest to the nth root of i."
return int(round(i ** (1./n)))
def relative_error(equation):
"Error between LHS and RHS of equation, relative to RHS."
(A, B, C, n) = equation
LHS = A ** n + B ** n
RHS = C ** n
return abs(LHS - RHS) / RHS
if __name__ == '__main__':
print(tests())
print("Searching beal(500, 100)")
print(beal(500, 100))
print("Finding Simpson-esque near-solutions to Fermat's Equation")
def s(b, p): print('{0}^{3} + {1}^{3} = {2}^{3}'.format(*simpsons(b, p)))
s(range(1000, 2000), [11, 12, 13])
s(range(3000, 5000), [12])
print("Searching beal_modp(500, 100)")
print(beal_modp(500, 100))