Blogdemaths - Théorème des deux carrés de Fermat

#########################################
# blogdemaths.wordpress.com
##########################################
# Théorème des deux carrés de Fermat
#
# https://blogdemaths.wordpress.com/2014/12/24/le-cadeau-de-noel-de-fermat/
#
# L'algorithme suivant renvoie la décomposition
# en somme de deux carrés d'un nombre premier p
# congru à 1 modulo 4
#
#
##########################################

def racine_carree_de_moins_un(p):
    """ renvoie un entier tel que u^2 = - 1 mod (p)"""

    # On commence par déterminer q et s tels que
    # p-1 = q * 2^s avec q impair

    s = p - 1
    e = 0
    while s % 2 == 0:
        s = s// 2
        e += 1

    #On détermine un entier n qui n'est pas un carré modulo p
    n=2
    while pow(n,(p-1)//2,p)==1: #Critère d'Euler
        n+=1

    # Le code suivant est issu de:
    # http://eli.thegreenplace.net/2009/03/07/computing-modular-square-roots-in-python/
    x= pow(-1,(s+1)//2,p)
    b = pow(-1,s,p)
    g = pow(n,s,p)
    r = e

    while True:
        m=0
        t=b

        for m in range(0,r):
            if t==1:
                break
            t= pow(t,2,p)

        if m==0:
            return x

        gs = pow(g,2**(r-m-1),p)
        g = (gs*gs) % p
        x = (x * gs) % p
        b = (b * g) % p
        r = m



def decomposition(x,y,m,p):
    """ A partir d'une décomposition x^2 + y^2 = m*p
       renvoie x1, y1, m1 tels que x1^2 + y_1^2 = m1*p
       avec m1 < m """

    r = x%m
    if r > m/2:
        r-=m

    s = y%m
    if s > m/2:
        s-=m

    m1 =  (r**2 + s**2)//m
    x1 = abs((r*x + s*y)//m)
    y1 = abs( (r*y - s*x)//m )

    return x1, y1, m1



p = int(input("Entrez un nombre premier p congru à 1 modulo 4:\n"))
u = racine_carree_de_moins_un(p)
m = (u**2 + 1)//p
#print("Etapes successives de l'algorithme:\n")
#print("{0}*{1} = {2}^2 + {3}^2".format(m,p,u,1))


x1,y1,m1 = decomposition(u,1,m,p)
#print("{0}*{1} = {2}^2 + {3}^2".format(m1,p,x1,y1))
while(m1>1):
        x1,y1,m1 = decomposition(x1,y1,m1,p)
        #print("{0}*{1} = {2}^2 + {3}^2".format(m1,p,x1,y1))


print("\nLa décomposition de {0} en une somme de deux carrés est:".format(p))
print("{0} = {1}^2 + {2}^2".format(p,x1,y1))