\\ Pollard Rho Faktorisierung
\\ N die zu faktorisierende Zahl,
\\ c die Konstante in der quadratischen Iterationsfunktion
\\ x der Startwert der Iteration,
\\ iter der Zweierlogarithmus der Anzahl der durchzuführenden
\\ Iterationen. Also: Erhöhung um 1 verdoppelt die Rechenzeit!
pollard_rho(N,c,x,iter)=
{
local(potenz1,potenz2,i,y);
potenz1=1; potenz2=2;
x=Mod(x,N);y=x;c=Mod(c,N);
while(iter-=1,
for(i=potenz1,potenz2,
y=y^2+c;
if(gcd(lift(x)-lift(y),N)>1,
print(iter);
return(gcd(lift(x)-lift(y),N));
);
);
x=y; potenz1*=2; potenz2*=2;
);
print("Iterationsgrenze überschritten!\n");
}
\\ Pollard (p-1)-Test
\\ N ist die zu faktorisierende Zahl,
\\ B Schranke: es werden alle Primzahlpotenzen, die unterhalb B liegen,
\\ miteinander multipliziert zu einer großen Zahl k
\\ Wenn p ein Primfaktor von N ist und jede Primzahlpotenz in der Primfaktorzerlegung
\\ von (p-1) unterhalb B liegt, dann ist offenbar (p-1) ein Teiler von B.
\\ Also wird x^k-1 durch p teilbar sein.
\\ Also wird x^k-1 einen gemeinsamen Teiler mit N besitzen
\\ Natürlich brauchen wir aber x^k-1 nur mod N auszurechnen.
\\ Beispiel: mit factor(nextprime(random(10^15))-1) erhielt ich die Zahl 884614672101733,
\\ deren größter Primfaktor 191 ist.
pollard_pminus1(N,B,x)=
{
local(p,q,k);
p=q=2;k=1;
while(p<B, while((q*=p)<B,print(q)); k*=(q/p); p=q=nextprime(p+1);print(p));
print(p);
print(k);
return(gcd(lift(Mod(x,N)^k)-1,N));
}
\\ findsmooth(n,B) sucht eine zufällige Primzahl in der Größenordnung n, für die die
\\ Faktorisierung von p-1 einen größten Primfaktor kleiner als B besitzt.
\\ findsmooth(10^15,1000) führt schnell zum Erfolg.
\\ findsmooth(10^20,10000) führt schnell zu p=1723846016192119513 mit 3527 als max. Faktor
\\ von p-1. Damit ist die 40stellige Zahl N=12507506288189026238869578837961025190149, die
\\ durch p teilbar ist, mit dem Aufruf pollard_pminus1(N,5000,2) faktorisierbar.
\\ findsmooth(10^25,100000) führt in 30sec zu p=3107286682619449216120777 mit 27967 als
\\ maximalem Primfaktor von p-1.
findsmooth(n,B) =
{
while(1,
m=factor(nextprime(random(n))-1);
v=matsize(m);
if(m[v[1],1]<B,
return(factorback(m)+1);
)
);
}
\\ Fermat Faktorisierung
\\ Zunächst zufällige Erzeugung von Relationen
\\ N ist die zu faktorisierende Zahl
\\ primebasesize ist die Größe der Primfaktorbasis
\\ iter ist die Anzahl der Versuche, zufällig glatte Quadrate bzgl. der Basis zu finden.
Fermat(N,primebasesize,iter)=
{
local(primenums,a,avec,b,bvec,b1,i,j,m,m2,primefactors,counter);
primenums=vector(primebasesize);
primenums[1]=2;
for(i=2,primebasesize,primenums[i]=nextprime(primenums[i-1]+1));
m=matrix(primebasesize,0);
avec=vectorv(0);
bvec=vectorv(0);
for(j=1,iter,
\\ für ein zufälliges a wird geschaut, ob b=a^2 glatt ist.
a=random(N); b1=b=lift(Mod(a,N)^2);
primefactors=vectorv(primebasesize);
for(i=1,primebasesize,
while((b1>=primenums[i]) && (!(b1%primenums[i])),
primefactors[i]+=1; b1 /= primenums[i];
);
);
if(b1==1,
avec=concat(avec,a);bvec=concat(bvec,b);m=concat(m,primefactors);
if((counter+=1)>=primebasesize*1.2,break);
);
);
\\ Die Matrix m enhält jetzt in jeder Spalte die Primzahlzerlegung eines gefundenen b,
\\ m2 enhält die Relationen
\\ avec und bvec enthalten die gefundenen a,b=a^2
m2=Mod(m,2);
m2=lift(matker(m2));
print("iter=",iter);
print("counter=",counter);
print(matsize(m2)[2]," Relationen gefunden.");
counter=matsize(m2)[1];
for(j=1,matsize(m2)[2],
a=lift( prod(i=1,counter,Mod(avec[i],N)^m2[i,j]));
b=lift(Mod(sqrtint(prod(i=1,counter,bvec[i]^m2[i,j])),N));
print(gcd(a+b,N));
);
\\ return([m,m2,avec,bvec]);
}
\\ Fermat(1555966091,30,20000), Fermat(610476462373,90,30000)
\\ Als naechstes folgt die Fermat-Faktorisierung mit Hilfe von Kettenbrüchen.
\\ Wesentlich dabei ist, daß die Quadrate b jedenfalls < Wurzel(N) sind, also
\\ nur halb so lang und deshalb mit einer deutlich kleineren Faktorbasis faktorisierbar.
\\ Wesentlich ist der Bhaskara Brouncker Algorithmus.
next_square_prime(n,N)=
{
local(p);
p=nextprime(n);
while(kronecker(N,p)!=1,p=nextprime(p+1));
return(p);
}
Fermat_CFrac(N,primebasesize,iter)=
{
local(primenums,a,avec,b,bvec,b1,i,j,m,m2,primefactors,counter);
primenums=vector(primebasesize);
primenums[1]=-1;primenums[2]=next_square_prime(1,N);
for(i=3,primebasesize,primenums[i]=next_square_prime(primenums[i-1]+1,N));
print(primenums);
m=matrix(primebasesize,0);
avec=vectorv(0);
bvec=vectorv(0);
A_CF=vector(3);B_CF=vector(3);C_CF=vector(3);P_CF=vector(3);i_CF=vector(3,i,i);
sqrt_CF=sqrtint(N);
A_CF[i_CF[2]]=sqrt_CF;
B_CF[i_CF[1]]=0;B_CF[i_CF[2]]=sqrt_CF;
C_CF[i_CF[1]]=1;C_CF[i_CF[2]]=N-sqrt_CF*sqrt_CF;
P_CF[i_CF[1]]=1;P_CF[i_CF[2]]=sqrt_CF;
vorzeichen=1;vorzeichenwechsel=-1;
for(j=1,iter,
\\ Der Bhaskara-Brouncker Algorithmus wirft jetzt a,b heraus mit a^2=b mod N,
\\ wobei die b dem Betrage nach < als 2sqrt(N) sind. Man sieht auch sofort, daß
\\ man nur solche Primzahlen in der Faktor-Basis braucht, für die (n/p)=1
A_CF[i_CF[3]]=(sqrt_CF+B_CF[i_CF[2]])\C_CF[i_CF[2]];
B_CF[i_CF[3]]=A_CF[i_CF[3]]*C_CF[i_CF[2]]-B_CF[i_CF[2]];
C_CF[i_CF[3]]=C_CF[i_CF[1]]+A_CF[i_CF[3]]*(B_CF[i_CF[2]]-B_CF[i_CF[3]]);
P_CF[i_CF[3]]=(P_CF[i_CF[1]]+A_CF[i_CF[3]]*P_CF[i_CF[2]])%N;
tmp=i_CF[1];i_CF[1]=i_CF[2];i_CF[2]=i_CF[3];i_CF[3]=tmp;vorzeichen*=vorzeichenwechsel;
a=P_CF[i_CF[1]]; b1=b=vorzeichen*C_CF[i_CF[1]];
\\ print(a," ",b);
primefactors=vectorv(primebasesize);
if(b<0,primefactors[1]+=1;b1=-b1);
for(i=2,primebasesize,
while((b1>=primenums[i]) && (!(b1%primenums[i])),
primefactors[i]+=1; b1 /= primenums[i];
);
);
if(b1==1,
avec=concat(avec,a);bvec=concat(bvec,b);m=concat(m,primefactors);
if((counter+=1)>=primebasesize*1.2,break);
);
);
\\ Die Matrix m enhält jetzt in jeder Spalte die Primzahlzerlegung eines gefundenen b,
\\ m2 enhält die Relationen
\\ avec und bvec enthalten die gefundenen a,b=a^2
m2=Mod(m,2);
m2=lift(matker(m2));
print("counter=",counter);
print(matsize(m2)[2]," Relationen gefunden.");
counter=matsize(m2)[1];
for(j=1,matsize(m2)[2],
a=lift( prod(i=1,counter,Mod(avec[i],N)^m2[i,j]));
b=lift(Mod(sqrtint(prod(i=1,counter,bvec[i]^m2[i,j])),N));
print(gcd(a+b,N));
);
\\ return([m,m2,avec,bvec]);
}