Problem 6 https://www.youtube.com/watch?v=Y30VF3cSIYQ

problem6.py

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)