Blatt 5, Aufgabe 3

Hier die Musterlösung:

def find_squares( p):
    """
    Return the solution (x,y) of 0<x<=y of p=x^2+y^2 if it exists,
    otherwise throw an exception. Input must be prime p = +1 mod 4.
    """
    assert p%4 == 1, 'Error: %d must be a prime = 1 mod 4' % p
    # Find a solution of x^2=-1 mod p
    # using Wilson's theorem
    x = 1
    for j in range(2,(p-1)/2):
        x *= j
        x = x%p

    # Modify x if necessary so that p > x > p/2
    if 2*x < p:
        x = p - x

    # compute r_l as in the preceding theorem
    a=p; b=x
    while b*b > p:
        r=a%b; a=b; b=r
          
    a = int( math.sqrt(p-b*b))
    return (a,b) if a < b else (b,a)

can = [p for p in primes(10^6, 10^6+100) if p%4 == 1]
can

# Lösungen berechnen
sol = [ (p,find_squares( p)) for p in can]

# Lösungen ausgeben und überprüfen
for p,(a,b) in sol:
    print p, a, b, p == a^2+b^2