Nous allons dans ce document présenter quelques résultats pour expliquer progressivement une procédure probabiliste de factorisation de grands nombres entiers, dite procédure de Pollard Floyd. 

La présente page fournit des explications incomplètes sur le fonctionnement de cette procédure, et se veut une invitation à panacher mathématique théorique et mathématique expérimentale, en utilisant les qualités propres aux deux types de logiciels que sont Geoplan-Geospace et Maple.
Les lecteurs sont donc conviés à une lecture active, de sorte à prolonger ce travail et à en combler les vides explicatifs.

Une figure Géoplan mettra en œuvre ce procédé de Pollard-Floyd en fin de document. Cette figure permet de décomposer des nombres composites de l'ordre de 10¹⁰.  Il va de soi que le présent document n'a donc pas la prétention de flirter avec les records actuels ; ce document vise en revanche à vulgariser des mathématiques qui agitent de plus en plus la mathématiciens et informaticiens, agitation qui accompagne l'invasion de notre quotidien par les technologies digitales.

A - Le "paradoxe" probabiliste des anniversaires :  

Une population de n individus se répartit en deux groupes:  le groupe A de ceux qui sont seuls à posséder leur date de naissance et le groupe B de ceux qui partagent leur date de naissance avec au moins une autre personne dans la population. Nous ne connaissons pas le cardinal de A et B. Nous allons donc aborder le problème d'un point de vue probabiliste, en nous basant sur le seul fait qu'il y a 365 (=N) jours dans une année (on oublie volontairement les 29 février et on considèrera que les individus de la population ont des dates de naissance indépendantes les unes des autres: pas de jumeaux par exemple).

Évaluation de la probabilité qu'un individu tiré au hasard dans la population soit né le même jour qu'un ou plusieurs autres individus déjà sondés.
Notons P(X=k) la probabilité que les k premières personnes sondées aient chacune leur anniversaire à des jours différents.
La première personne sondée peut avoir son anniversaire n'importe quel jour.
La seconde personne sondée a 364/365 ou 1-1/365 chances d'avoir son anniversaire un autre jour, la troisième 363/365 ou 1-2/365, et ainsi de suite...
Donc P(X=1)=1, P(X=2)= 1-1/365 et plus généralement, pour tout entier k>1, on a la relation de récurrence 
P(X=k)=P(X=k-1).(1-(k-1)/365) puisque P(X=k)=(1-1/365).(1-2/365).(1-3/365)...(1-(k-1)/365)  .
Il va de soi que P(X=366)=0, et a fortiori pour tout k >366, P(X=k)=0, ce dont rend bien compte la formule exhibée.

Donc la probabilité de l'événement contraire, à savoir de l'événement "au moins deux personnes  parmi les k sondées ont même date d'anniversaire" est P( X=k ) = 1- P(X=k) . Cette suite est évidemment strictement croissante entre 1 et 365, puis constante égale à 1. Un item de Géoplan (indice du premier terme nul d'une suite) permet de répondre au problème d'antécédent, à savoir calcule k en fonction d'une probabilité p.
Plus précisément, la figure suivante permet de répondre aux questions du type suivant, où p est choisi  arbitrairement entre 0 et 1: "à partir de combien de personnes interrogées a-t-on une probabilité supérieure à p d'avoir au moins un doublon?"  Les réponses surprennent (d'où le recours à l'expression paradoxe des anniversaires).  
Par exemple pour p=1/2, on trouve k = 23:  il suffit de sonder 23 personnes pour avoir plus d'une chance sur deux qu'une date d'anniversaire au moins soit partagée.
Plus généralement, la figure suivante permet d'explorer une généralisation du paradoxe des anniversaires de N=365 à N quelconque, (on parle alors des N attributs possibles de la population) et de p=1/2 à p probabilité quelconque.

Pour piloter le nombre N d'attributs : touche N, puis flèches du clavier ; 
pour piloter p: touche P, puis flèches du clavier, ou directement à la souris

De plus, par développement de Taylor à l'ordre 1,    est approximativement égal à  .

On peut donc approcher P(X=k), qui en toute rigueur est égal à (1-1/365)(1-2/365)(1-3/365)...(1-(k-1)/365), par

p' =    

et donc  l'événement "au moins deux personnes sondées parmi les k ont même date d'anniversaire" a une probabilité 

qui peut être approchée par p =  .

Cette dernière expression fournit, c'est la moindre des choses, une suite croissante en la variable entière k>1. Inversement, pour p probabilité fixé, il est aisé de montrer que k doit être à peu près égal à   ; en particulier pour p=1/2, on a k qui est de l'ordre de  , fonction de N qui est asymptotiquement équivalent à  , et donc de l'ordre de  .

B - le calcul d'une longueur de cycle suivant une idée de Floyd

Soit f une fonction d'un ensemble E de cardinal fini dans lui-même.
On définit grâce à cette fonction f et un élément u₀ de E une suite u récurrente d'ordre 1 : u(0)=u₀, et  u(n+1)=f(u(n)) ; 
La suite u est bien évidemment périodique à partir d'un certain rang, disons le rang q. En effet, la suite u peut être considérée comme une fonction de N dans E, et pour des raisons de cardinal, u ne peut pas être injective. Soit donc (i,j) un couple d'indices tels que u_i=u_j, et i<j. Puisque la suite est récurrente d'ordre 1, cela suffit à établir l'existence d'une partie périodique dont la longueur p divise nécessairement j-i , la période commençant au plus tard à partir du rang i  (autrement dit i>q-1). Si u est périodique à partir du rang q, on dira que q est la longueur de la pré-période de la suite u. La proposition suivante fonde un algorithme efficace de calcul à la fois d'un majorant de q, un multiple de la longueur p de la période, et aussi et surtout pour la section suivante d'un couple (i,j), i<j, avec u_i=u_j.

La démonstration  de cette première proposition peut être introduite par une image, de sorte à soutenir l'intuition: le parcours d'une suite périodique à partir d'un certain rang peut être dessiné comme la lettre ρ : une partie apériodique ou plus précisément pré-périodique, à savoir la queue de la lettre ρ, suivie d'une partie périodique, la boucle de la lettre  ρ.

Démonstration Soit donc une suite périodique à partir d'un rang q. Les deux suites u et v peuvent donc être interprétées comme deux coureurs cyclistes,- un amateur qui va deux fois moins vite qu'un pro- partant en même temps à une distance de q unités d'un vélodrome et se retrouvant tôt ou tard au même point à l'intérieur du vélodrome. Si on admet dans un premier temps l'existence d'un premier rang i de rencontre, alors il est clair que le pro aura parcouru un nombre entier de périodes de plus que l'amateur, puisque les deux coureurs ont d'abord chacun parcouru à leur rythme une seule fois la pré-période, et, juste avant leur rencontre, la même portion du vélodrome. L'existence d'un premier rang i de rencontre est due au fait qu'une fois le coureur lent engagé dans la période, il a un certain nombre d'unités de longueur d' "avance" sur l'autre, avance qu'il se fait grignoter unité après unité. (L'existence d'un rang de rencontre ne serait pas acquise si la "suite rapide" w allait 3 fois plus vite que la suite lente, c'est-à-dire si on avait posé w(n)=(fofof)(w(n-1))=u(3n)  et non w(n)=(fof)(w(n-1))=u(2n) ): cela ferait deux unités de longueur grignotées par itération, et il faudrait alors discuter selon les parités des longueurs q de la pré-période et p de la période.
En résumé: si u est périodique à partir d'un rang q, alors il existe un plus petit rang de rencontre i>q-1 pour lequel u(i)=w(i), et i est nécessairement multiple de la longueur p de la période de la suite u. 
Par suite w(i+p)=((f ²)^p)(w(i)) = ((f ²)^p)(u(i))=((f ^p)²)(u(i))=(f ^p)(u(i+p))=(f ^p)u(i)=u(i+p), et par récurrence w(i+np)=u(i+np) pour tout n entier.
On vient de montrer que E(f) contient tous les multiples de p à partir de i. 
On rappelle que i est défini comme le plus petit  (ou, dans notre version coureur cycliste, le premier) rang de rencontre, et non comme le plus petit entier multiple de p qui soit supérieur au sens large à q. 
Définissons donc  i' comme le plus petit multiple de p qui soit supérieur au sens large à q; puisqu'il a été établi que i est lui même multiple de p, i' est par définition inférieur ou égal à i. Si on arrive à démontrer que u(i')=w(i'), alors i sera également par définition inférieur ou égal à i'. Supposons u(i') différent de w(i') ; on aurait alors, exactement comme plus haut, en remplaçant les i par des i' et les signes égal à par des signes différent de, w(i'+np) différent de u(i'+np) pour tout n entier, en particulier, comme pour un certain n, i'+np=i, on aurait w(i) différent de u(i), ce qui est contradictoire avec l'existence établie de i et sa caractérisation, contradiction qui prouve que u(i') = w(i') et donc que i'=i. Finalement, l'ensemble des multiples de p qui sont supérieurs à q est inclus dans E(f). On laisse au lecteur démontrer l'inclusion inverse, en s'inspirant de ce qui a été dit en début de démonstration sur l'indice i.

C - Application: procédure probabiliste de factorisation de grands nombres entiers, dit procédé de Pollard-Floyd

On considère pour un entier n l'ensemble F={0,1,...., n-1}muni de sa structure canonique d'anneau et g une fonction de F dans F, qui pourrait par exemple être définie par 
g(x) = x²+1 modulo n pour tout x de F.
La fonction g va tenir simultanément deux rôles:
    1 -  g, couplée à  un premier élément  v(0) de F aléatoirement choisi, va définir de manière pseudo-aléatoire la suite des éléments de F à sonder. Plus précisément, le premier élément sondé sera v(0), le deuxième sera v(1) = g(v(0)), etc...
    2 -  de plus, g tient un rôle d'attribution, c'est-à-dire card(g(F)) va tenir le rôle du nombre 365 dans le paradoxe des anniversaires; 
l'ensemble A regroupe les éléments x de F qui ne partagent  leur attribut g(x) avec aucun autre élément de F; une condition nécessaire pour que x soit élément de A est que la suite v initialisée par v(0)=x ait une pré-période de longueur q différente de 1;
l'ensemble B regroupe les éléments y de F qui partagent leur attribut g(y) avec au moins un autre élément de F; une condition suffisante pour que y soit élément de B est que la suite v initialisée par v(0)=y ait une pré-période de longueur égale à 1 (q=1); autrement dit, il suffit que y soit l'unique élément de la pré-période de la suite v qu'il initialise pour que y soit élément de B; en effet, la longueur p de la période étant toujours strictement positive, donc ici nécessairement supérieure ou égale à q, l'élément y'=v(p) est dans la période , donc y' est distinct de y; de plus, y et y' ont même attribut: g(y)=v(1)=v(p+1)=g(v(p))=g(y').

D'après l'analyse faite plus haut du paradoxe des anniversaires, et sous réserve que le double rôle de g n'entame pas trop le bien-fondé de l'hypothèse "le sondage peut être considéré comme aléatoire", on commence à mettre les chances de notre côté pour obtenir un doublon après avoir sondé un nombre d'éléments de l'ordre de  . qui est bien sûr inférieur à   (égalité si g est bijective, ce qui ne serait pas avantageux), c'est-à-dire on peut espérer trouver un couple (i,j) avec i<j=O( ) et  v(i)=v(j).
Ceci étant dit, soit à présent N un (grand) nombre entier composite dont on désire exhiber un facteur non trivial. Dans la suite, N admettra comme facteur non trivial l'entier n et N=nK. 
Remarquons tout de suite que la fonction g  ne peut pas convenir pour un algorithme de factorisation, pour la bonne raison que n est a priori inconnu, au même titre que tout facteur non trivial de N.
On peut en revanche tenter notre chance avec f définie par f(x)=x²+1 modulo N.  Cela contraint à considérer non plus l'ensemble F={0,1,...., n-1}, mais l'ensemble E={0 , ..., N-1}. C'est donc cette fois-ci la fonction  f  et non plus g qui tient simultanément deux rôles:  le rôle de co-définition ( via la relation de récurrence u(n+1)=f(u(n)) ) de la marche du sondeur, une fois un élément u(0) choisi aléatoirement dans E  et, de plus, le rôle d'attribution (card(f(E)) est à notre problème ce que 365 est au paradoxe des anniversaires;  la partition de E en A U B vaut également, où A est désormais l'ensemble des éléments x de E qui ne partagent  leur attribut f(x) avec aucun autre élément de E; une condition nécessaire pour que x soit élément de A est que la suite u initialisée par u(0)=x ait une pré-période de longueur q différente de 1;
l'ensemble B regroupe les éléments y de E qui partagent leur attribut f(y) avec au moins un autre élément de E; une condition suffisante pour que y soit élément de B est que la suite u initialisée par u(0)=y ait une pré-période de longueur égale à 1 (q=1);
En utilisant la factorisation supposée N=nK, il est aisé de démontrer que pour tout x de E que g(x)=f(x) modulo n ; en particulier, on a,  pour tout i>0, v(i) = u(i)  modulo n si l'on prend u(0) = v(0).
Puisque l'on peut espérer, moyennant une procédure de sondage, trouver un couple (i,j) avec i<j=O( ) et  v(i)=v(j), on peut également espérer trouver un couple (i,j) avec i<j=O( ) et u(i)=u(j)  modulo n , autrement dit un couple qui permette d'écrire u(i)-u(j)=nL  pour un certain entier L.  On a alors pgcd(u(i)-u(j),N)=pgcd(nK, nL)=n*pgcd(K,L). 
Un tel couple (i,j) suscite donc la discussion suivante:  
si L est non nul , alors pgcd(K,L) sera compris entre 1 et min(K,L);
      si pgcd(K,L)<min(K,L), alors pgcd(u(i)-u(j),N) est un facteur non trivial de N que l'on exhibera; 
      sinon, il est toujours envisageable de changer l'expression de la fonction f ou le premier élément u(0) à sonder.
sinon, L=0 et alors pgcd(u(i)-u(j),N)=pgcd(0,N)=N qui est facteur trivial de N.
Enfin, si contrairement aux attentes pgcd(u(i)-u(j),N)=1, c'est que le raisonnement fait de puis le départ avec n facteur non trivial vaut aussi avec n=1...et pour les nombres N qui sont premiers.
Remarques:
- On a vu plus haut que l'on majorait   par   . Si la fonction g concernée était bijective, on aurait égalité   , ce qui ne serait pas vraiment avantageux; cependant, à l'autre extrême, si card(g(F))=1, c'est-à-dire si g est constante, on obtiendrait tout de suite des couples (i,j) ad hoc, mais pour conclure trivialement de manière bien vaine que N divise N (cas L=0 dans la discussion menée plus haut); pour rendre efficace la méthode, il faudra chercher une fonction g non bijective (ce qui proscrit g(x)=ax+b avec a et n premiers entre eux) mais tout de même plus riche qu'une vaine fonction constante.
- les facteurs d'un entier N marchant par paires, on est tenté d'avancer que O( ) peut être remplacé par O( ). Toutefois, bien que la borne   habituellement recommandée dans la littérature consacrée à cette méthode fournisse des résultats plus que satisfaisants, nous restons dubitatifs quant à l'argument. Lorsque l'on teste la procédure sur des nombres tels que :
                         N=648 637 207 ( dont les deux seuls facteurs non triviaux sont 7 207 et 90 001), alors c'est au bout de 2848 itérations que la factorisation aboutit avec la fonction f(x)=x²+1 modulo N ; or    =25 468 environ, et  =160 environ. Pire, si  N=5 777 076 197 (dont les deux seuls facteurs non triviaux sont 71 789 et 80 473), alors c'est au bout de 36 763 itérations que la factorisation aboutit avec la fonction f(x)=3x²+1 modulo N ; or  = 76 007 environ, et =275 environ. Ces deux exemples montrent que la borne en   a la vertu d'accélérer la vitesse de réponse, mais empêche parfois de faire aboutir certaines factorisations délicates. Les fonctions explicitées ci-dessus ont co-définies la marche de sondage avec le premier élément sondé u(0)= 1.

Une fois clairement établi le but , à savoir trouver un couple (i,j) tel que i<j et u(i)=u(j), reste à déterminer quel algorithme permet d'y arriver efficacement. C'est là que l'astuce de Floyd nous est précieuse: nous sommes dans les conditions d'application de la propriété de la section B.  Plutôt que de tester un à un une liste exhaustive de tous  les couples (i,j) de E², plutôt que de garder en mémoire les u(i),  il est préférable d'extraire la suite w définie dans la propriété de la section B pour (accessoirement) déterminer un multiple de la longueur p du cycle et un majorant de la longueur q de la pré-période et (essentiellement) pour repérer un couple (i,j) ad hoc (avec donc j=2i), couple qui pourra nous fournir un facteur de N, sauf si pgcd(u(i)-u(j),N)=N ou 1.
Jugez  par vous-même de la rapidité de la réponse d'éventuelle factorisation, même pour des nombres entiers de l'ordre 10¹⁰; la figure ci-dessous fournit donc avec certitude de nombreuses décompositions et, en cas d'échec de décomposition, rend fortement suspect de primalité le nombre considéré. On ne peut annoncer avec certitude que le nombre est alors premier, mais après vérification avec un logiciel de calcul formel (voir annexe), les tests menés avec la seule fonction f₁(x)=x²+1 fournissent des nombres premiers de l'ordre du million en ne se trompant "presque jamais". 
Pour doubler l'attaque, on se donne une seconde fonction polynomiale, choisi expérimentalement, à savoir  f₂(x)=3x²+1. 
On laisse au lecteur la liberté de redéfinir la double attaque: un double clic sur l'une des deux expressions
 f1(x)=  ou f2(x)= ,  une redéfinition des fonctions dans la fenêtre "EcritMath" ouverte, et, pour en sortir, un simple appui sur "appliquer" puis sur le bouton .
Il n'est pas interdit de tripler l'attaque , étant donné la rapidité de la réponse fournie par l'algorithme.

Pour piloter l'entier N à factoriser:  flèches du clavier. 

Il faudrait encore examiner avec un oeil critique l'hypothèse selon laquelle nos fonctions d'attaque, qui ont le double rôle d'attribution et de co-définition de la suite des éléments sondés, sont "suffisamment aléatoires", voire, pourquoi pas, "mieux qu'aléatoires". 
La fonction f₁(x)=x²+1 revient systématiquement dans la littérature consacrée à la méthode de factorisation de Pollard Floyd. Cependant des nombres composites modestes comme 125 = 5³ ou 169=13² résistent à l'attaque de f₁. On démontre en théorie algébrique des nombres que les seuls diviseurs premiers possibles de x²+1 pour x entier sont 2 (mais nous ne considérons ici que les nombres impairs) et des nombres premiers du type 1+4k;  5 et 13 rentrent dans cette catégorie; est-ce une coïncidence fortuite ?
L'étude générale des liens entre une fonction polynomiale f et l'ensemble des entiers composites qui lui résistent est un problème difficile.
Une recherche moins ambitieuse et néanmoins difficile consiste à chercher des fonctions f puissantes, c'est-à-dire des fonctions qui échouent rarement à exhiber une décomposition de N composite, à l'instar de cette fonction f₁(x)=x²+1. R. Guy (cf. références) a commencé à étudier ce problème, et a pour cela introduit la notion d'epact d'un nombre premier n relativement à une fonction f définie sur Z/nZ, qui est le plus grand entier k pour lequel il existe un entier s avec card{s,f(s),...,f^k(s)}=k+1. Autrement dit, l'epact est à une unité près la taille maximale de la somme p+q (somme des longueurs de la pré-période et de la période), lorsque f est fixé et que l'on joue sur le paramètre s d'initialisation de la suite. Les choix optimaux de fonctions d'attaque sont d'après R. Guy ceux qui minimisent cette somme p+q. Il apparaît alors (pour une seconde raison) que les fonctions affines sont à proscrire, tout comme certaines fonctions quadratiques, telles f(x)=x²-2.
R. Guy conjectura que la somme p+q maximale de f(x)=x²+1 a une croissance de l'ordre de O( ).

Une autre conjecture à ce jour non prouvée est que il y a une infinité de nombres premiers de la forme x²+1.
Remarquons qu'à l'exception peut-être des publications de R. Guy référencées ci-dessous et que nous n'avons pas eu loisir de lire, nous n'avons pas eu à préciser que les facteurs n de N devaient être premiers. 

Signbalons enfin que les fondateurs de l'informatique (Von Neumann en tête)  avaient repéré le problème des pré-périodes et périodes, en poursuivant un autre but. Ils cherchaient à mettre au point des générateurs pseudo-aléatoires de nombres, et pour ces générateurs pseudo-aléatoires, (où u(0) est qualifié de "graine" - seed en anglais), il est préférable que la pré-période soit minimale.

Références:
Une approche très riche mi- mathématique, mi- informatique de l’arithmétique des grands nombres : 
Prime Numbers, A computational Perspective, de Richard Crandall et Carl Pomerance.

De plus, ce livre fournit des références sur la méthode de factorisation de Pollard-Floyd :
J. Pollard: A Monte Carlo Method for factorization.Nordisk Tidskr. Informationsbehandling (BIT) ,15:331-334, 1975.
R. Guy :  How to factor a number.  Proceedings of the fifth Manitoba Conference on Numerical Mathematics ; vol 16 of Congress Congressus Numerantium, pp. 49-89, 1976
R. Brent et J . Pollard : Factorization of the eighth Fermat Number. Math. Comp., 36 : 627-630 , 1981.

_________________________________________________________________________________________

Annexe: On peut vérifier que les nombres suspects de primalité sont effectivement premiers ou non avec un logiciel de calcul formel, comme par exemple avec Maple:

with(numtheory):
isprime(7633123); 
nextprime(7633123);

true
7633147

isprime(7633139); factorset(7633139);

false
{71,107509}

La procédure de Pollard-Floyd peut être testée grâce à d'autres procédures, comme celles, implémentées dans le logiciel Maple, numtheory[divisors], ou encore  isprime, dont on peut afficher l'algorithme comme suit:
  interface(verboseproc=3):print(isprime);

proc (n) local btor, nr, p, r; option remember, system, 
`Copyright (c) 1993 Gaston Gonnet, Wissenschaftliches Rechnen, ETH Zurich. All rights reserved.`;

if not type(n,integer) then 
      if type(n,numeric) then ERROR(`argument must be an integer`) 
      else RETURN('isprime(n)') 
      fi
 fi; 
if n < 2 then false
elif has(`isprime/w`,n) then true
elif igcd(2305567963945518424753102147331756070,n) <> 1 then false
elif n < 10201 then true
elif igcd(84969694892334181105323399091873499659260625866489327366
115454263422038932707693909090694773095091375097869171186680288
614993338250976823867229837379629630667576741311267365789364407
881571869698937306331130664786204486249492573240226273954373636
390387526081667586612559568346306972204475122988482222285500626
837863425199602259963013159456444700647206966217504772445289159
27867113,n) <> 1 then false
elif n < 1018081 then true
 else nr := igcd(408410100000,n-1); 
nr := igcd(nr^5,n-1); 
r := iquo(n-1,nr); 
btor := modp(power(2,r),n);
if `isprime/cyclotest`(n,btor,2,r) = false or irem(nr,3) = 0 and `isprime/cyclotest`(n,btor,3,r) = false or irem(nr,5) = 0 and `isprime/cyclotest`(n,btor,5,r) = false or irem(nr,7) = 0 and `isprime/cyclotest`(n,btor,7,r) = false then RETURN(false) 
fi;
if isqrt(n)^2 = n then RETURN(false) fi; 
for p from 3 while numtheory[jacobi](p^2-4,n) <> -1 do od; evalb(`isprime/TraceModQF`(p,n+1,n) = [2, p]) 
fi
end;