Crible quadratique (factorisation)

Soyez le premier à donner votre avis sur cette source.

Vue 13 786 fois - Téléchargée 479 fois

Description

Mon code factorise des grands nombres assez rapidement (60 chiffres->20 minutes à 1.7GHz) en utilisant le MPQS qui est dérivé du crible quadratique.
D'habitude, le java est utilisé pour la programmation orientée objet mais je l'ai utilisé pour ces routines rapides de grands nombres.

Source / Exemple :


import java.math.BigInteger;
import java.util.BitSet;
class  MPQS{
	  /*
    Racine carrée entière.
    Renvoie le nombre maximal dont le carré est inférieur ou égal à n
    [Méthode dite de Newton]

  • /
static BigInteger sqrt(BigInteger n) { int nb = n.bitLength() ; BigInteger a ; a = BigInteger.valueOf(1); a = a.shiftLeft((nb+1)/2) ; for (int i=0;i<=1001;i++){ a = ((n.divide(a)).add(a)).shiftRight(1); } return a ; } /* Calcul du symbole de legendre, noté bizarrement (n/p) Pour p premier on a : - (n/p) = 0, ssi p divise n. Et sinon, - (n/p) = 1 ,ssi n est un résidu quadratique modulo p. - (n/p) = -1 ,ssi n n'est pas un résidu quadratique modulo p.
  • /
private static int legendre(BigInteger n, int p) { int a = n.remainder(BigInteger.valueOf(p)).intValue() ; int m = p ; int t = 1 ; while (a != 0) { int k = m % 8 ; boolean zyva = k == 3 || k == 5 ; while (a % 2 == 0) { a /= 2 ; if (zyva) t = -t ; } int swap = a ; a = m ; m = swap ; int x = a % 4, y = m % 4 ; if (x == 3 && x == y) t = -t ; a = a % m ; } if (m == 1) return t ; else return 0 ; } private static int multMod(int x, int y, int m) { long w = (long)x * (long)y ; return (int)(w % m) ; } static int powMod(int x, int k, int m) { if (k == 0) return 1 ; else { int r = powMod(x, k/2, m) ; r = multMod(r, r, m) ; if (k % 2 == 0) return r ; else return multMod(x, r, m) ; } } /* Calculer une solution de l'équation x^2 = n (mod p). Important: On suppose que cette solution existe, c'est dire que n est un résidu quadratique modulo p. [Algorithme de Shanks-Tonelli.]
  • /
private static int sqrtMod(BigInteger n, int p) { int a = n.remainder(BigInteger.valueOf(p)).intValue() ; if (p % 4 == 3) { return powMod(a, (p+1)/4, p) ; } else { int d ; do { d = 2+Integer.parseInt(""+Math.round(Math.random()*(p-2))) ; } while (legendre(BigInteger.valueOf(d), p)==1) ; int s = 0 ; int t = p-1 ; while (t % 2 == 0) { s++ ; t >>=1; } int A = powMod(a, t, p) ; int D = powMod(d, t, p) ; int m = 0 ; int ADpowM = A ; int twoPowI = 1 ; int D2 = D ; for (int exp = 1 << s-1 ; exp !=0 ; exp >>>= 1, twoPowI <<= 1, D2 = multMod(D2, D2, p)) { if (powMod(ADpowM, exp, p) == p-1) { m += twoPowI ; ADpowM = multMod(ADpowM, D2, p) ; } } return multMod(powMod(a, (t+1)/2, p), powMod(D, m/2, p), p) ; } } static int legendre(BigInteger a,BigInteger p){ BigInteger ret=a.modPow(p.shiftRight(1),p); if (ret.equals(p.subtract(BigI(1))))return -1; //System.out.println(ret); return ret.intValue(); } private static BigInteger sqrtMod(BigInteger n, BigInteger p) { BigInteger a = n.remainder(p); if (p.remainder(BigI(4)).equals(BigI(3))) { return a.modPow((p.add(BigI(1))).shiftRight(2), p) ; } else { int d ; do { d = 2+Integer.parseInt(""+Math.round(Math.random()*1000)) ; } while (legendre(BigI(d),p)==1) ; int s = 0 ; BigInteger t = p.subtract(BigI(1)); while (t.mod(BigI(2)).equals(BigI(0))) { s++ ; t=t.shiftRight(1); } BigInteger A =a.modPow(t,p) ; BigInteger D =BigI(d).modPow(t,p) ; int m = 0 ; BigInteger ADpowM = A ; int twoPowI = 1 ; BigInteger D2 = D ; for (int exp = 1 << s-1 ; exp !=0 ; exp >>>= 1, twoPowI <<= 1, D2 =D2.multiply(D2).remainder(p)) { if (ADpowM.modPow(BigI(exp),p).equals(p.subtract(BigI(1)))) { m += twoPowI ; ADpowM =ADpowM.multiply(D2).remainder(p); } } return (a.modPow(t.add(BigI(1)).shiftRight(1), p).multiply(D.modPow(BigI(m>>1), p))).remainder(p) ; } } static byte log2(int n){ // renvoie le log en base 2 int k=n; byte log=1; while (k!=1) { k>>=1; // k=k/2 log++; } return log; } static boolean prima(int i){ if (i==2)return true; if (i==3)return true; if (i==5)return true; final int rac=(int)Math.sqrt(i)+1; int k; for (k=2;k!=rac;k++) if (i%k==0)break; if (k==rac)return true; return false; } static BigInteger BigI(int i){ return(BigInteger.valueOf(i)); } static BigInteger f(long i,BigInteger racine,BigInteger n){ // calcule (racine+i)^2-n BigInteger k=racine.add(BigInteger.valueOf(i)); return k.multiply(k).subtract(n); } static BigInteger Pow(int x,int k){ //renvoie x^k if (k == 0) return BigI(1); else { BigInteger r=Pow(x, k/2) ; r=r.multiply(r); if (k % 2 == 0) return r ; else return r.multiply(BigI(x)) ; } } public static BigInteger nextprime(BigInteger a){ final BigInteger deux=BigI(2); if (a.equals(deux)) return BigI(3); BigInteger essai=a.add(deux); if (a.remainder(deux).equals(BigI(0))) essai=essai.subtract(BigI(1)); do{essai=essai.add(deux);}while (!essai.isProbablePrime(80)); return essai; } public static void main(String args[]){ long deb=System.currentTimeMillis(); final int taille=300000; //final int base=41903; final double ln2=Math.log(2); final int sup=30; final int min=10; final int tailleM1=taille-1; final BigInteger un=BigI(1); final BigInteger zero=BigI(0); final int multiplier=1; int tailprem=1,j,j2,k,nblisse=-1,nb=0,nbfaux=0,premi,racPrem,tmp; BigInteger n=new BigInteger("676292275716558246502605230897191366469551764092181362779759"), // entier à factoriser inter,ent,racine,fact,premj,carre,a,q,b,facto,ainv,b0,b1,h; n=n.multiply(BigI(multiplier)); final double lnN=Math.log(n.doubleValue()); final double temp=4*Math.exp(Math.sqrt(lnN*Math.log(lnN)/12)); final int base=(int)(temp*Math.log(temp)); System.out.println(base); byte expo,ajout,borne; q=sqrt(sqrt(BigI(2).multiply(n)).divide(BigI(taille))); racine=sqrt(n).add(un); System.out.println(n.divide(BigI(multiplier))); for (int i=2;i<=base;i++) if (prima(i)&&legendre(n,i)==1)tailprem++; int[] prem=new int [tailprem+1]; tailprem=1; prem[0]=-1; for (int i=2;i<=base;i++) if (prima(i)&&legendre(n,i)==1) { prem[tailprem]=i; tailprem++; } final int li=tailprem+sup,col=tailprem; final int liM1=li-1; BigInteger[] tab =new BigInteger[li]; BigInteger[] tabQ=new BigInteger[li]; byte[][] matrexp=new byte [li][col]; BitSet[]matr=new BitSet[li]; { int[] tsqrt=new int [col]; byte[] log2=new byte[col]; byte[] w=new byte [taille]; byte[] zeros=new byte [taille]; for (int i=0;i<taille;i++)zeros[i]=0; for (int i=0;i<li;i++)matr[i]=new BitSet(col); for (int i=min;i!=tailprem;i++){ premi=prem[i]; tsqrt[i]=sqrtMod(n,premi); log2[i]=(byte)Math.ceil(Math.log(premi)/ln2); } System.out.println(tailprem+" "+base); while (nblisse<li){ // tant que l'on a pas assez de résidus nb++; q=nextprime(q); while (legendre(n,q)!=1)q=nextprime(q); a=q.pow(2); b0=sqrtMod(n,q); h=(n.subtract(b0.pow(2))).divide(q); b1=(h.multiply((BigI(2).multiply(b0)).modInverse(q))).mod(q); b=(b0.add(b1.multiply(q))).remainder(a); if ((b.remainder(BigI(2))).equals(BigI(0)))b=a.subtract(b); borne=(byte)(racine.multiply(BigI(taille)).bitLength()-log2(base)); System.arraycopy(zeros,0,w,0,taille); for(k=min;k!=tailprem-1;k++){ premi=prem[k]; ajout=log2[k]; ainv=a.modInverse(BigI(premi)); j =((ainv.multiply(BigI( tsqrt[k]).subtract(b))).mod(BigI(premi))).intValue(); j2=((ainv.multiply(BigI(-tsqrt[k]).subtract(b))).mod(BigI(premi))).intValue(); if (j<j2){ tmp=j2; j2=j; j=tmp; } j2=j-j2; while (j<taille){ w[j]+=ajout; w[j-j2]+=ajout; j+=premi; } } for (int i=0;i<taille;i++){ for(;w[i]<borne;i++)if (i==tailleM1)break; if (i==taille)break; if (w[i]>=borne) { // si le nb à l'air lisse(tout ses facteurs sont plus petits que b nblisse++; facto=(a.multiply(BigI(i))).add(b); ent=((facto.pow(2)).subtract(n)).divide(a); // on le factorise matr[nblisse].clear(); if (ent.compareTo(zero)<0){ ent=ent.negate(); matrexp[nblisse][0]=1; matr[nblisse].set(0); }else{ matrexp[nblisse][0]=0; } tab [nblisse]=facto; tabQ[nblisse]=q; j=1;expo=0;premj=BigI(2); while(true){ if (ent.remainder(premj).equals(zero)){ ent=ent.divide(premj); expo++; if (ent.equals(un)){ matrexp[nblisse][j]=(byte)expo; // Il faut écrire le dernier exposant if ((expo&1)==1)matr[nblisse].set(j); System.arraycopy(zeros,0,matrexp[nblisse],j+1,tailprem-j-1); for (j++;j<col;j++)matr[nblisse].clear(j); break; // La factorisation est complète } } else{ matrexp[nblisse][j]=(byte)expo; if ((expo&1)==1)matr[nblisse].set(j); expo=0; j++; if (j==tailprem){ nblisse--; nbfaux++; break; // Les nombres premiers de la base ont été parcouru, le nombre n'est pas lisse } premj=BigInteger.valueOf(prem[j]); } } if (nblisse==(li-1)){ // On a trouvé toutes les relations suffisantes nblisse++; break; } } } if (nb%20==0) System.out.println("Fini : "+nblisse+"("+nbfaux+")="+(nblisse*100)/(tailprem+sup)+"% "); } } System.out.println(""); System.out.println("J'ai fini le crible avec "+nbfaux+" qui sont faux"); System.out.println("Criblage fini en "+(System.currentTimeMillis()-deb)+" ms"); System.out.println("Un nombre lisse toutes les "+((long)nb*(long)taille)/li+" valeurs"); System.out.println("Crible sur "+nb*2+" polynomes"); BitSet[]T=new BitSet[li]; BitSet bonPiv=new BitSet(li); j=0; for (int i=0;i<li;i++){ T[i]=new BitSet(li); bonPiv.set(i); for (int l=0;l<li;l++) T[i].clear(l); } for (int i=0;i<li;i++) T[i].set(i); // T devient une matrice identité int piv,rech; for (int m=0;m<col;m++){ // on fait le pivot de Gauss piv=0; while(true){ try{ while (!matr[piv].get(m)) piv++; } catch(Exception e){piv=-1;}; if (piv==-1) break; if (bonPiv.get(piv))break; piv++; }; if (piv==-1)continue; bonPiv.clear(piv); j2=-1; while (true){ if (j2==liM1)break; j2++; try{ while (!matr[j2].get(m))j2++; } catch(Exception e){j2=-1;}// On a dépassé les bornes du tableau if (j2==-1)break; if (j2!=piv){ matr[j2].xor(matr[piv]); T[j2].xor(T[piv]); } } } for(int i=0;i!=li;i++)matr[i]=new BitSet(0); for(int i=0;i!=col;i++)T [i]=new BitSet(0); System.out.println("Pivot de Gauss termine"); int []ligneExp=new int [col]; for (int i=col;i<li;i++){ BigInteger prodQ=un; ent=BigI(1);inter=BigI(1); for (int l=0;l<col;l++)ligneExp[l]=0; for (int l=0;l<li;l++) if (T[i].get(l)){ ent=(ent.multiply(tab[l])).mod(n); for (int m=0;m<col;m++)ligneExp[m]+=matrexp[l][m]; prodQ=(prodQ.multiply(tabQ[l])).remainder(n); } for (int l=0;l<col;l++) inter=(inter.multiply(Pow(prem[l],ligneExp[l]/2))).remainder(n); fact=n.gcd(ent.subtract(inter.multiply(prodQ))); if((!fact.equals(un))&&(!fact.equals(n))){ if (fact.equals(BigI(multiplier)))continue; if ((n.divide(fact)).equals(BigI(multiplier)))continue; if (fact.mod(BigI(multiplier)).equals(zero)) System.out.println((fact.divide(BigI(multiplier)))+"*"+n.divide(fact)); else System.out.println((fact)+"*"+((n.divide(fact)).divide(BigI(multiplier)))); break; } } System.out.println((System.currentTimeMillis()-deb)+" millisecondes"); } }

Conclusion :


Pour l'explication, regarder ici : http://www.enseignement.polytechnique.fr/profs/informatique/Luc.Maranget/IF/qs/sujet.html et http://mersennewiki.org/index.php/MPQS.

Pour l'implémentation, j'ai dû "découper en tranches" tout mon crible : au lieu de faire un tableau de 5 000 000 de case, j'ai préféré
(et je suis même obligé) de cribler par tranches de 100 000.
J'ai eu un problème : comment trouver la case qu'il faut éliminer?
Réponse : j=prem[k]-(nb*taille-tabj[k])%prem[k]
prem étant le tableau de nb premiers, tabj est le tableau des 1ère racines carrées.
taille=100000 qui est la taille d'une tranche
nb est le nombre de tranches

Cette méthode à un avantage : on n'a pas à dire "Je vais cribler jusqu'à x" mais "Je crible tant que je n'ai pas assez de résidus".
Car je ne sais pas comment calculer x.

Pour le MPQS, j'ai mis une formule pour la taille de la base.
Vous pouvez modifier la taille de l'intervalle. Cela peut accélérer le programme.

Codes Sources

A voir également

Ajouter un commentaire Commentaires
salut ,bah je vous remercie pour ce programme ,j'aimerais bien si vous pouvez m'offrir la documentation expliquant en detail qs ou mpqs pour voir deja ou se situe la difference clairement et je vous serai tres reconnaissant si vous avez la meme implementation pour qs ou mpqs en C meme pour les petits entier j vasi le transformer ensuite avc la biblio gmp pour les grands nombre ,merci.
Messages postés
1
Date d'inscription
dimanche 19 novembre 2006
Statut
Membre
Dernière intervention
27 juillet 2008

peu robuste et lent
Messages postés
20
Date d'inscription
mardi 11 octobre 2005
Statut
Membre
Dernière intervention
13 mars 2007

Ce qui me dégoute, c'est cette page : http://www.enseignement.polytechnique.fr/profs/informatique/Luc.Maranget/IF/qs/main002.html
J'arrive en 10s à factoriser T20 : 5 fois plus lent.
J'arrive en 12s à factoriser T30 : 3 fois plus lent.
J'arrive en 3min30 à factoriser F7 : 7 fois plus lent.
J'arrive en 4min40 à factoriser T40 : 9 fois plus lent.
J'arrive en 11min à factoriser T45 : 5 fois plus lent.
Et je ne vois pas comment optimiser encore plus.
Messages postés
20
Date d'inscription
mardi 11 octobre 2005
Statut
Membre
Dernière intervention
13 mars 2007

Oui mais je n'ai jamais vu l'intérêt de cette méthode. Ca rapporte à la même procedure, faire un objet BigInteger dedans revient à le mettre en paramètre. Donc inutile de le faire. Autant ne pas se compliquer la vie.
Messages postés
2835
Date d'inscription
lundi 11 août 2003
Statut
Modérateur
Dernière intervention
16 décembre 2006
10
Effectivement c'est tout sauf du Java ton code, m'enfin utile pour ceux qui voudront utiliser cette méthode ... ;-)

D'autant qu'il ne t'aurais pas fallu grand chose (30 sec) pour le mettre en objet ... ;-)

++

Vous n'êtes pas encore membre ?

inscrivez-vous, c'est gratuit et ça prend moins d'une minute !

Les membres obtiennent plus de réponses que les utilisateurs anonymes.

Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.

Le fait d'être membre vous permet d'avoir des options supplémentaires.