Skip to content

Instantly share code, notes, and snippets.

@skaslev
Last active January 9, 2020 09:20
Show Gist options
  • Save skaslev/52b59474f8b20fb3ff654a87ddbf1181 to your computer and use it in GitHub Desktop.
Save skaslev/52b59474f8b20fb3ff654a87ddbf1181 to your computer and use it in GitHub Desktop.
from sys import stdout
from primefac import primefac
from math import sqrt
# -- "any sufficiently well-commented lisp program contains an ML program in its comments"
# def solution (p a b : N) := a^2 + b^2 = (ab + 1) p^2
# def sym (p a b : N) : solution p a b -> solution p b a
# def sol0 (p : N) : solution p 0 p
# def generator (p a b : N) (p > 0) (a <= b) : solution p a b -> solution p b (b * p^2 - a)
# Thm: If a, b such that a^2 + b^2 = (ab + 1) p^2 then
# b^2 + (bp^2-a)^2 = (b(bp^2-a) + 1) p^2
# e.g. a' = b and b' = bp^2-a are also solution of the equation
# Proof.
# Simplifying the left hand side:
# b^2 + (bp^2-a)^2 = a^2 + b^2 - 2abp^2 + b^2 p^4 = b^2 p^4 + (1 - ab) p^2
# and the rhs reduces to:
# (b(bp^2-a) + 1) p^2 = b^2 p^4 + (1 - ab) p^2
# Qed.
# Q{0} = 0
# Q{1} = p
# Q{n+2} = p^2 Q{n+1} - Q{n}
# --
# Q{2} = p^3
# Q{3} = p^5 - p
# Q{4} = p^7 - 2p^3
# Q{5} = p^9 - 3p^5 + p
# Q{6} = p^11 - 4p^7 + 3p^3
# Q{7} = p^13 - 5p^9 + 6p^5 - p
# Q{8} = p^15 - 6p^11 + 10p^7 - 4p^3
# Q{9} = p^17 - 7p^13 + 15p^9 - 10p^5 + p
#
# (Q{n+1}) (0 1) (Q{n} )
# ( ) = ( ) ( )
# (Q{n+2}) (-1 p^2) (Q{n+1})
def Q(n, p):
if n == 0: return 0
if n == 1: return p
return p**2*Q(n-1,p)-Q(n-2,p)
def genp(N, p):
a, b = 0, p
while a < N and b < N:
print a, b, sqrt(1.0*(a**2+b**2)/(a*b+1))
a, b = b, b * p**2 - a
if b == 0:
break
def gen(N, P):
for p in range(P+1):
genp(N, p)
def log(x, base):
try:
return math.log(x, base)
except:
return float('nan')
def search(N):
for a in range(N):
print '\r', a,
stdout.flush()
for b in range(a,N):
if (a**2 + b**2) % (a * b + 1) != 0:
continue
# For all b, (0, b) is a solution since
# (0^2+b^2)/(0b+1) = b^2
if a == 0: # and b == p:
continue
# For all a, (a, a^3) is a solution since
# (a^2+a^6)/(a^4+1) = a^2(1+a^4)/(1+a^4) = a^2
if a**3 == b: # e.g. a == p and b == p**3:
continue
p = int(sqrt((a**2 + b**2) / (a * b + 1)))
# For all p, (p^3, p^5-p) is a solution since
# (p^6+(p^5-p)^2)/(p^8-p^4+1) = p^2(p^8-p^4+1)/(p^8-p^4+1) = p^2
if a == p**3 and b == p**5-p:
continue
# For all p, (p^5-p, p^7-2p^3) is a solution since
# p^2(p^8-2p^4+1+p^12-4p^8+4p^4)/(p^12-2p^8-p^8+2p^4+1) = p^2(p^12-3p^8+2p^4+1)/(p^12-3p^8+2p^4+1)= p^2
if a == p**5-p and b == p**7-2*p**3:
continue
# Keep finding new solutions by iterating
# a, b = b, b * p**2 - a
if (a == p**7-2*p**3 and b == p**9-3*p**5+p or
a == p**9-3*p**5+p and b == p**11-4*p**7+3*p**3 or
a == p**11-4*p**7+3*p**3 and b == p**13 - 5*p**9 + 6*p**5 - p or
a == p**13 - 5*p**9 + 6*p**5 - p and b == p**15 - 6*p**11 + 10*p**7 - 4*p**3):
continue
print '\r',
print(a, b, (a**2 + b**2) / (a * b + 1), sqrt((a**2 + b**2) / (a * b + 1)), log(b, a))
print(list(primefac(a)))
print(list(primefac(b)))
print("")
if __name__ == '__main__':
gen(10**9, 20)
# search(10**6)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment