IMPLEMENTATION MODULE SUPFACT; (* Transportabilite: Le plus grand nombre premier du fichier 'SUPFACT.PRI' *) (* doit etre inferieur a 2^(Beta/2), ou Beta est le nombre de bits *) (* utilises pour representer un entier. BornePremExp2 doit etre inferieur *) (* a P*P, ou P est le plus petit nombre premier superieur au plus grand *) (* nombre premier du fichier 'SUPFACT.PRI'. *) (* Dans cette version, Beta = 32 et BornePremExp2 = 2^32. *) FROM storage IMPORT ALLOCATE; FROM MYINOUT IMPORT WriteString, WriteLn, ReadLongInt, ReadCard, SetInput, GetInput, WriteCard, EOF, Read, ReadLn, CloseInput, OpenAsInput; (* Les fonctions GetInput, CloseInput, OpenAsInput de MYINOUT remplace les *) (* fonctions Open, Close precedemment utilisees sur le systeme SUN *) FROM FileIO IMPORT FILE; FROM ErrorSystem IMPORT Result; FROM BasicIO IMPORT Ok; FROM SUPINT IMPORT SuperInteger, Create, Delete, Affect; FROM SUPCONV IMPORT IntToSup, SupToCard; FROM SUPARITH IMPORT Even, Negative, Zero, Smaller, Greater, Equal, Add, Subtract, Multiply, Divide, DivideSupInt, Quotient, Modulo, PowerSupInt, GCD, MultMod, PowerMod; FROM SUPRAND IMPORT SupRandom; FROM CHRONO IMPORT Chrono, TimeUnit, TimeFormat, InitChrono, ValChrono; (* FROM Heure IMPORT Heure; *) (* car CHRONO l'utilise et Heure est ecrit en C *) CONST LenSupInt = 4; (* Taille d'un SuperInteger mod N. *) NEssaisC = 100; (* Nombre d'essais pour tenter de *) (* prouver qu'un nombre est compose *) (* pendant une factorisation. *) NomFchNbPrem = "/home/matawin/stochas/sentiers/mcs/supfact.pri"; (* Nom du fichier des nombres premiers *) (* inferieurs a 2^16. *) VAR F : Factors; (* Facteur courant dans la liste L. *) FirstFact : BOOLEAN; (* =TRUE si le 1er fact. n'a pas *) (* encore ete trouve. *) Timer : Chrono; (* Constantes. *) Sup0, Sup1, Sup2, BornePremExp2, (* Borne sup. sur le plus grand nombre *) (* premier du fichier au carre (=2^32).*) (* Variables vraiment globales. *) N, (* Cofacteur restant a decomposer. *) NTemp, ProduitFact, (* Variables locales a certaines procedures, mais declarees globales *) (* pour sauver des appels a Create. *) b, m, AMoins1, NbPremExp2, Reste, Germe, T1, Temp, TestTemp, ProdTemp : SuperInteger; c : CHAR; (****************************************************************************) PROCEDURE IsProbablyPrime ( A : SuperInteger; k : CARDINAL (* Nombre de tentatives pour prouver que *) (* A est compose. *) ) : BOOLEAN; (* Retourne FALSE si A n'est pas premier (A est alors compose). *) (* Retourne TRUE si l'algorithme n'a pu determiner que A est compose. *) (* Si le resultat est TRUE, la probabilite que A soit un nombre compose *) (* et qu'il n'ait pas ete detecte comme tel est < 2 ** (- k). *) (* Voir H.S. Wilf, "Algorithms and Complexity", p. 158 et p. 164. *) VAR i, q : LONGINT; essai : CARDINAL; Compose : BOOLEAN; BEGIN IF Even ( A ) THEN IF Equal (A,Sup2) THEN RETURN TRUE ELSE RETURN FALSE END; END; (* A est impair. *) Subtract ( A, Sup1, AMoins1 ); Affect ( AMoins1, m ); (* m = A - 1 est pair. *) q := 0; REPEAT DivideSupInt ( m, 2, m, i ); INC ( q ) UNTIL NOT Even ( m ); (* AMoins1 = m * 2^q. *) essai := 0; REPEAT SupRandom ( Sup2, AMoins1, b ); (* 2 =< b =< A - 1 *) PowerMod ( b, m, A, T1 ) ; (* T1 = b**m mod A *) Compose := NOT ( Equal ( T1, Sup1 ) OR Equal ( T1, AMoins1 ) ); (* b^m <> 1, b^m <> A-1 ( mod A ). *) i := 1; WHILE (i <= q - 1) AND Compose DO MultMod ( T1, T1, A, T1 ); (* T1 = b^(m*2^i) mod A *) Compose := NOT Equal ( T1, AMoins1 ); (* b^(m*2^i) <> A-1 (mod A). *) INC ( i ) END; INC ( essai ) UNTIL (essai = k) OR Compose; RETURN NOT Compose END IsProbablyPrime; (****************************************************************************) PROCEDURE IsPrime ( A : SuperInteger; k : CARDINAL ) : BOOLEAN; (* Idem a IsProbablyPrime, mais verifie d'abord si A est divisible par *) (* un nombre premier inferieur a 2^16. *) VAR NbPrem : CARDINAL; a, NbPremExp2 : LONGCARD; Reste : INTEGER; CurrFchIn : FILE; BEGIN IF Even ( A ) THEN IF Equal (A,Sup2) THEN RETURN TRUE ELSE RETURN FALSE END; END; (* A est impair. *) GetInput (CurrFchIn); OpenAsInput (NomFchNbPrem); IF NOT (Result() = Ok) THEN WriteString ("** ERREUR IsPrime ** "); WriteString ("Fichier des nombres premiers 'supfact.pri' introuvable."); WriteLn; HALT END; ReadCard(NbPrem); IF Smaller ( A, BornePremExp2 ) THEN (* A < 2^32, on peut prouver sa primalite ou sa non-primalite. *) SupToCard ( A, a ); NbPremExp2 := 1; WHILE (NOT EOF()) AND ( a > NbPremExp2 ) DO ReadLn; IF (a MOD NbPrem) = 0 THEN CloseInput; SetInput (CurrFchIn); RETURN FALSE END; (* A MOD NbPrem <> 0. *) NbPremExp2 := NbPrem * NbPrem; Read(c); IF NOT EOF() THEN ReadCard(NbPrem); END; END; (* A est premier. *) CloseInput; SetInput (CurrFchIn); RETURN TRUE ELSE (* A >= 2^32. *) WHILE NOT EOF () DO ReadLn; DivideSupInt ( A, VAL(INTEGER,NbPrem), T1, Reste ); IF Reste = 0 THEN CloseInput; SetInput (CurrFchIn); RETURN FALSE END; (* A MOD NbPrem <> 0. *) Read(c); IF NOT EOF() THEN ReadCard(NbPrem); END; END END; CloseInput; SetInput (CurrFchIn); (* A n'est divisible par aucun des nombres premiers inf. a 2^16. *) RETURN IsProbablyPrime ( A, k ) END IsPrime; (****************************************************************************) PROCEDURE NewFact ( A : SuperInteger; VAR L : Factors ); (* Place un nouveau facteur A de N a la fin de la liste L. *) (* NOTE : NTemp doit etre egal a N / A. *) BEGIN IF FirstFact THEN FirstFact := FALSE; NEW ( L ); F := L ELSE NEW ( F^.Next ); F := F^.Next END; WITH F^ DO Create ( Factor, LenSupInt ); Affect ( A, Factor ); Multiplicity := 0; Prime := TRUE; REPEAT INC ( Multiplicity ); Affect ( NTemp, N ); Divide ( N, A, NTemp, Reste ); UNTIL NOT Zero ( Reste ) OR Equal ( N, Sup1 ) END; END NewFact; (****************************************************************************) PROCEDURE NewCoFact ( A : SuperInteger; VAR L : Factors ); (* Place le cofacteur compose A a la fin de la liste L. *) BEGIN IF FirstFact THEN FirstFact := FALSE; NEW ( L ); F := L ELSE NEW ( F^.Next ); F := F^.Next END; Divide ( N, A, NTemp, Reste ); WITH F^ DO Create ( Factor, LenSupInt ); Affect ( A, Factor ); Multiplicity := 0; Prime := FALSE; REPEAT INC ( Multiplicity ); Affect ( NTemp, N ); Divide ( N, A, NTemp, Reste ) UNTIL NOT Zero ( Reste ) OR Equal ( N, Sup1 ); END END NewCoFact; (****************************************************************************) PROCEDURE TrialDivision ( VAR L : Factors ); (* Tente de diviser N par chacun des nombres premiers d'un fichier, *) (* jusqu'a 65536. Place les facteurs trouves dans L. A la fin, N est *) (* soit egal a 1, auquel cas il a ete completement factorise, soit *) (* egal au cofacteur restant a decomposer. Ce cofacteur est premier *) (* s'il est inferieur a 65536 * 65536. *) VAR NbPrem : CARDINAL; Reste : LONGINT; CurrFchIn : FILE; BEGIN GetInput (CurrFchIn); OpenAsInput (NomFchNbPrem); IF NOT (Result() = Ok) THEN WriteString ("** ERREUR TrialDivision ** "); WriteString ("Fichier des nombres premiers 'supfact.pri' introuvable."); WriteLn; HALT END; IntToSup ( 1, NbPremExp2 ); ReadCard(NbPrem); WHILE (NOT EOF()) AND Greater ( N, NbPremExp2 ) DO ReadLn; DivideSupInt ( N, NbPrem, NTemp, Reste ); IntToSup ( NbPrem, NbPremExp2 ); IF Reste = 0 THEN (* NbPremExp2 = NbPrem est un nouveau facteur. *) NewFact ( NbPremExp2, L ); END; Multiply ( NbPremExp2, NbPremExp2, NbPremExp2 ); (* NbPremExp2 = NbPrem*NbPrem. *) Read(c); IF NOT EOF() THEN ReadCard(NbPrem); END; END; CloseInput; SetInput (CurrFchIn); IF Greater ( N, Sup1 ) AND Smaller ( N, NbPremExp2 ) THEN (* N est facteur premier du N original. *) IntToSup ( 1, NTemp ); NewFact ( N, L ); (* N = 1. *) END END TrialDivision; (****************************************************************************) PROCEDURE PollardMCF ( M : SuperInteger; T : LONGREAL; (* Temps max. de CPU en minutes. *) VAR L : Factors ); (* Tente de factoriser un entier M par la methode de Pollard *) (* (algorithme de type Monte-Carlo). Voir P.L. Montgomery, *) (* "Speeding the Pollard and Elliptic Curve Methods of Factorization", *) (* Math. of Computation, vol. 48, no 177, janvier 1987, pp. 243-264. *) (* Place les facteurs trouves a la fin de la liste L. Assume que tous *) (* les facteurs inf. a SQRT(BornePremExp2) ont ete enleves. *) CONST cmax = 2; (* Nb max. de tours avant d'appeler *) (* Check pour prendre un GCD. *) VAR e, f, i, j : INTEGER; d, (* Procedure Check. *) t1, t2, w, x, y, z, a1, a2, a3, a4 : SuperInteger; c : [0..cmax]; B : ARRAY [1..cmax] OF SuperInteger; Echec, FactComplete : BOOLEAN; PROCEDURE G ( x, y : SuperInteger ); (* Effectue l'affectation y := (x + 1) mod M. *) BEGIN Add ( x, Sup1, y ); IF NOT Smaller ( y, M ) THEN Subtract ( y, M, y ) END END G; PROCEDURE Check; (* Calcule le cofacteur et fait le backtracking quand on trouve un *) (* diviseur de M. *) VAR i, k : CARDINAL; dTemp : SuperInteger; (* Variable pour l'appel recursif a PollardMCF.*) BEGIN IF c > 0 THEN Affect ( Sup1, ProdTemp); FOR i := 1 TO c DO MultMod ( ProdTemp, B [i], M, ProdTemp ) END; GCD ( ProdTemp, M, Temp ); IF Greater ( Temp, Sup1 ) THEN REPEAT Affect ( M, d ); k := 0; FOR i := 1 TO c DO GCD ( B [i], M, Temp ); (* Temp = GCD ( B [i], M ). *) IF Greater ( Temp, Sup1 ) THEN IF NOT Equal ( M, d ) THEN GCD ( B [i], d, Temp ) END; (* Temp = GCD ( B [i], d ). *) IF Greater ( Temp, Sup1 ) THEN Affect ( Temp, d ) END; INC ( k ); Affect ( B [i], B [k] ); END; END; IF k > 0 THEN IF Equal ( M, d ) THEN Echec := TRUE ELSE (* M <> d. *) IF Greater ( BornePremExp2, d ) THEN (* d est premier. *) Quotient ( M, d, NTemp ); NewFact ( d, L ) ELSE IF IsProbablyPrime ( d, NEssaisC ) THEN (* d est probablement premier. *) Quotient ( M, d, NTemp ); NewFact ( d, L ) ELSE (* d est compose. *) Create ( dTemp, LenSupInt ); Affect ( d, dTemp ); PollardMCF ( dTemp, T, L ); Delete ( dTemp ) END END; c := k; END; END UNTIL k = 0; IF NOT Echec THEN IF Greater ( BornePremExp2, M ) THEN (* M est premier. *) IF M # N THEN Quotient ( N, M, NTemp ) ELSE IntToSup ( 1, NTemp ) END; NewFact ( M, L ); FactComplete := TRUE; ELSE IF IsProbablyPrime ( M, NEssaisC ) THEN (* M est probablement premier. *) IF M # N THEN Quotient ( N, M, NTemp ) ELSE IntToSup ( 1, NTemp ) END; NewFact ( M, L ); FactComplete := TRUE; ELSE (* Le cofacteur M est compose. *) IF ValChrono ( Timer, min ) >= T THEN (* Temps CPU depasse. *) NewCoFact ( M, L ); FactComplete := TRUE; RETURN END; (* Reduire les resultats intermediaires modulo M. *) Modulo ( t1, M, t1 ); Modulo ( t2, M, t2 ); Modulo ( w, M, w ); Modulo ( x, M, x ); Modulo ( y, M, y ); Modulo ( z, M, z ); Modulo ( a1, M, a1 ); Modulo ( a2, M, a2 ); Modulo ( a3, M, a3 ); Modulo ( a4, M, a4 ) END END END END END; IF (ValChrono ( Timer, min ) >= T) AND NOT FactComplete THEN (* Temps CPU depasse. *) NewCoFact ( M, L ); FactComplete := TRUE; RETURN END; c := 0; END Check; PROCEDURE InitPollardMCF; BEGIN Subtract ( M, Sup1, Temp ); (* Temp = M - 1. *) SupRandom ( Sup0, Temp, Germe ); Multiply ( Germe, Germe, z ); G ( z, z ); MultMod ( z, z, M, y ); c := 0; e := 1; f := 1; i := 1; Affect ( z, t1 ); Affect ( Germe, t2 ); G ( y, x ); END InitPollardMCF; BEGIN (* PollardMCF *) IF ValChrono ( Timer, min ) >= T THEN (* Temps CPU depasse. *) NewCoFact ( M, L ); RETURN END; Create ( w, LenSupInt ); Create ( x, LenSupInt ); Create ( y, LenSupInt ); Create ( z, LenSupInt ); Create ( t1, LenSupInt ); Create ( t2, LenSupInt ); Create ( a1, LenSupInt ); Create ( a2, LenSupInt ); Create ( a3, LenSupInt ); Create ( a4, LenSupInt ); Create ( d, LenSupInt ); FOR j := 1 TO cmax DO Create ( B [j], LenSupInt ) END; Echec := FALSE; FactComplete := FALSE; REPEAT InitPollardMCF; REPEAT Add ( t1, t2, a1 ); IF NOT Smaller ( a1, M ) THEN Subtract ( a1, M, a1 ) END; Add ( a1, z, a2 ); IF NOT Smaller ( a2, M ) THEN Subtract ( a2, M, a2 ) END; Multiply ( a1, z, a3 ); Multiply ( t1, t2, Temp ); Add ( a3, Temp, a3 ); Modulo ( a3, M, a3 ); Add ( y, a3, a4 ); MultMod ( a1, a4, M, a4 ); j := e + f; e := f; f := j; Affect ( t1, t2 ); Affect ( z, t1 ); REPEAT Subtract ( x, a2, w ); IF Negative ( w ) THEN Add ( w, M, w ) END; MultMod ( x, x, M, y ); Add (y, a3, TestTemp ); MultMod ( TestTemp, w, M, TestTemp ); Add ( TestTemp, a4, TestTemp ); IF NOT Smaller ( TestTemp, M ) THEN Subtract ( TestTemp, M, TestTemp ) END; (* TestTemp = (y + a3) * w + a4 *) (* = (x - t1)(x - t2)(x - z) (mod M). *) (* Test ( TestTemp ) *); INC ( c ); Affect ( TestTemp, B [c] ); IF (c = cmax) OR (ValChrono ( Timer, min ) >= T) THEN Check END; G ( y, x ); INC ( i ) UNTIL (i = f) OR FactComplete OR Echec; IF NOT ( FactComplete OR Echec ) THEN Add ( w, a2, z ); IF NOT Smaller ( z, M ) THEN Subtract ( z, M, z ) END END; UNTIL FactComplete OR Echec; UNTIL FactComplete; Delete ( t1 ); Delete ( t2 ); Delete ( w ); Delete ( x ); Delete ( y ); Delete ( z ); Delete ( a1 ); Delete ( a2 ); Delete ( a3 ); Delete ( a4 ); FOR j := 1 TO cmax DO Delete ( B [j] ) END END PollardMCF; (****************************************************************************) PROCEDURE Factorize ( A : SuperInteger; T : LONGREAL; (* Temps de CPU max. en minutes. *) VAR L : Factors (* Liste des facteurs trouves. *) ); (* Tente de trouver les facteurs premiers de A en appliquant successi- *) (* vement les procedures TrialDivision et PollardMCF. *) BEGIN InitChrono ( Timer ); FirstFact := TRUE; Affect ( A, N ); TrialDivision ( L ); IF Greater ( N , Sup1 ) THEN (* La factorisation n'est pas complete. *) IF IsProbablyPrime ( N, NEssaisC ) THEN (* N est facteur premier de A. *) IntToSup ( 1, NTemp ); NewFact ( N, L ) (* N = 1. *) ELSE (* N est compose. *) PollardMCF ( N, T, L ) END END; IF FirstFact THEN IntToSup ( 1, NTemp ); NewFact ( N, L ) END; F^.Next := NIL; (* Verifier le produit des facteurs obtenus. *) IntToSup ( 1, ProduitFact ); F := L; REPEAT WITH F^ DO PowerSupInt ( Factor, Multiplicity, Temp ); Multiply ( Temp, ProduitFact, ProduitFact ); F := Next END UNTIL F = NIL; IF NOT Equal ( A, ProduitFact ) THEN WriteString ('** ERREUR Factorize ** Produit des facteurs trouves <> A.'); WriteLn; HALT END END Factorize; (****************************************************************************) BEGIN (* Initialisation du module SUPFACT. *) Create ( Sup0, 0 ); IntToSup ( 0, Sup0 ); Create ( Sup1, 0 ); IntToSup ( 1, Sup1 ); Create ( Sup2, 0 ); IntToSup ( 2, Sup2 ); Create ( BornePremExp2, 2 ); PowerSupInt ( Sup2, 32, BornePremExp2 ); Create ( b, LenSupInt ); Create ( m, LenSupInt ); Create ( AMoins1, LenSupInt ); Create ( T1, LenSupInt ); Create ( ProduitFact, LenSupInt ); Create ( N, LenSupInt ); Create ( NTemp, LenSupInt ); Create ( NbPremExp2, 2 ); Create ( Germe, LenSupInt ); Create ( Reste, LenSupInt ); Create ( Temp, LenSupInt ); Create ( TestTemp, LenSupInt ); Create ( ProdTemp, LenSupInt ); END SUPFACT.