Blogdemaths

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

Dec 24th, 2014
652
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #########################################
  2. # blogdemaths.wordpress.com
  3. ##########################################
  4. # Théorème des deux carrés de Fermat
  5. #
  6. # https://blogdemaths.wordpress.com/2014/12/24/le-cadeau-de-noel-de-fermat/
  7. #
  8. # L'algorithme suivant renvoie la décomposition
  9. # en somme de deux carrés d'un nombre premier p
  10. # congru à 1 modulo 4
  11. #
  12. #
  13. ##########################################
  14.  
  15. def racine_carree_de_moins_un(p):
  16.     """ renvoie un entier tel que u^2 = - 1 mod (p)"""
  17.  
  18.     # On commence par déterminer q et s tels que
  19.     # p-1 = q * 2^s avec q impair
  20.    
  21.     s = p - 1
  22.     e = 0
  23.     while s % 2 == 0:
  24.         s = s// 2
  25.         e += 1
  26.    
  27.     #On détermine un entier n qui n'est pas un carré modulo p
  28.     n=2
  29.     while pow(n,(p-1)//2,p)==1: #Critère d'Euler
  30.         n+=1
  31.        
  32.     # Le code suivant est issu de:
  33.     # http://eli.thegreenplace.net/2009/03/07/computing-modular-square-roots-in-python/
  34.     x= pow(-1,(s+1)//2,p)
  35.     b = pow(-1,s,p)
  36.     g = pow(n,s,p)
  37.     r = e
  38.    
  39.     while True:
  40.         m=0
  41.         t=b
  42.  
  43.         for m in range(0,r):
  44.             if t==1:
  45.                 break
  46.             t= pow(t,2,p)
  47.  
  48.         if m==0:
  49.             return x
  50.  
  51.         gs = pow(g,2**(r-m-1),p)
  52.         g = (gs*gs) % p
  53.         x = (x * gs) % p
  54.         b = (b * g) % p
  55.         r = m
  56.        
  57.        
  58.  
  59. def decomposition(x,y,m,p):
  60.     """ A partir d'une décomposition x^2 + y^2 = m*p
  61.       renvoie x1, y1, m1 tels que x1^2 + y_1^2 = m1*p
  62.       avec m1 < m """
  63.  
  64.     r = x%m
  65.     if r > m/2:
  66.         r-=m
  67.  
  68.     s = y%m
  69.     if s > m/2:
  70.         s-=m
  71.    
  72.     m1 =  (r**2 + s**2)//m
  73.     x1 = abs((r*x + s*y)//m)
  74.     y1 = abs( (r*y - s*x)//m )
  75.  
  76.     return x1, y1, m1
  77.  
  78.  
  79.  
  80. p = int(input("Entrez un nombre premier p congru à 1 modulo 4:\n"))
  81. u = racine_carree_de_moins_un(p)
  82. m = (u**2 + 1)//p
  83. #print("Etapes successives de l'algorithme:\n")
  84. #print("{0}*{1} = {2}^2 + {3}^2".format(m,p,u,1))
  85.  
  86.  
  87. x1,y1,m1 = decomposition(u,1,m,p)
  88. #print("{0}*{1} = {2}^2 + {3}^2".format(m1,p,x1,y1))
  89. while(m1>1):
  90.         x1,y1,m1 = decomposition(x1,y1,m1,p)
  91.         #print("{0}*{1} = {2}^2 + {3}^2".format(m1,p,x1,y1))
  92.  
  93.        
  94. print("\nLa décomposition de {0} en une somme de deux carrés est:".format(p))
  95. print("{0} = {1}^2 + {2}^2".format(p,x1,y1))
RAW Paste Data