Google+ Followers

Translate

Post in evidenza

SimpleGNFS in PARI GP - parte 4

Ecco oggi ho un pò di tempo per mostrarvi praticamente una semplice versione di GNFS in PARI GP. Il simpleGNFS di Rosario Turco e niente pi...

sabato 4 febbraio 2017

SimpleGNFS in PARI GP - parte 4


Ecco oggi ho un pò di tempo per mostrarvi praticamente una semplice versione di GNFS in PARI GP. Il simpleGNFS di Rosario Turco e niente più teoria o chiacchiere. Solo applicazione.

Il programma simpleGNFS è costituito dalle seguenti parti:
- GNFS.gp    (programma principale)
- Comb.gp    (routine utility per combinazioni semplici usato anche in CFRAC)
- SUMMOD2.gp (routine di algebra lineare)
- SetupGNFS.txt (setup del GNFS)

Tutti i file devono essere posizionati nella stessa directory.

SetupGNFS.txt, Comb.gp e SUMMOD2.gp sono inclusi in GNFS.txt col comando \r \, dove il nome della directory va modificata in base a dove si posizionano i file.

Comb.gp lo abbiamo usato anche in CFRAC e qui ci aiuta per due tipi di lista di combinazioni semplici:
- quella che ci porta ai primi ideali A={(ri,pi)}
- quella che ci porta ad una soluzione della matrice X; nel seguito E è la matrice già modulo 2 rispetto a X

SUMMOD2.gp trova, invece, le combinazioni lineari  dipendenti che portano ad una soluzione non banale nel gcd.

Il file SetupGNFS.txt è da configurare prima di un lancio e per ogni n che volete fattorizzare. Tutti i parametri presenti in tale file vanno configurati. I valori presenti inizialmente sono utili per fattorizzare ad esempio n=45113.

L'algoritmo è didattico e andrebbe ottimizzato. Ad esempio Comb.gp ha la funzionalità CombList() che difficilmente ce la fa ad elencare tutte le combinazioni di 50 elementi (per un fatto di memoria). POtrete fare qualche allocatemem() per aumentare lo stack ma fino ad un certo punto. Inoltre è consigliato un linguaggio come il C, per avere una maggiore velocità. Infine è da ricordare maggiore è n da fattorizzare più risorse servono (memoria e tempo). GNFS per numeri veramente grandi può anche metterci diversi mesi a fattorizzare.

Si ipotizza,poi, che il polinomio f(x) venga scelto a mano e scritto nel SetupGNFS.txt

Vedremo, in qualche altro post, che è possibile, invece, automatizzare anche la scelta del polinomio f(x), che qui viene inserito nel SETUP.

Una volta posizionato in una directory tutti i file e modificato in GNFS.gp il nome della directory in \r,  il lancio di simpleGNFS avviene in tre passi:
1) scegliere n. Ad esempio
n=45113
2) \r /GNFS.gp
3) simpleGNFS(n)

Il programma simpleGNFS(45113), con la configurazione di prima, inizia a fare calcoli e propone anche un fattore algebrico A suggerendo di usare solo 20 elementi (volteAtoR è settato a 2 cioè A è 2 volte R). Maggiore è tale valore volteAtoR più lento è il programma e maggiore è la memoria necessaria. Si parte col valore 2 e si aumenta nel setup solo se non si riesce a fattorizzare.

Tipicamente con A={(ri,pi)}} si sceglieranno una ventina di coppie (ri,pi) con pi distinti e in qualche caso duplicato se questo serve per completare il numero di elementi.

Ad esempio suggerisco di inserire per A la seguente lista (fate copia e incolla):
[[0,2],[6,7],[13,17],[11,23],[26,29],[18,31],[19,41],[13,43],[1,53],[46,61],[2,67],[6,67],[44,67],[50,73],[23,79],[47,79],[73,79],[28,89],[62,89],[73,89]];

Ecco l'output del programma

gp > simpleGNFS(45113)

Numero composto

Inizio SETUP

Polinomio configurato in SETUP

factor base R di 10 elementi

== R ==
factor base R
List([2, 3, 5, 7, 11, 13, 17, 19, 23, 29])

factor base algebrico A. Esempio:
[0,2],[2,2],[4,2],[6,2],[8,2],[10,2],[12,2],[14,2],[16,2],[18,2],[20,2],[22,2],[24,2],[26,2],[28,2],[30,2],[32,2],[34,2],[36,2],[38,2],[40,2],[42,2],[44,2],[46,2],[48,2],[50,2],[52,2],[54,2],[56,2],[58,2],[60,2],[62,2],[64,2],[66,2],[68,2],[70,2],[72,2],[74,2],[76,2],[78,2],[80,2],[82,2],[84,2],[86,2],[88,2],[90,2],[6,7],[13,7],[20,7],[27,7],[34,7],[41,7],[48,7],[55,7],[62,7],[69,7],[76,7],[83,7],[90,7],[13,17],[30,17],[47,17],[64,17],[81,17],[11,23],[34,23],[57,23],[80,23],[26,29],[55,29],[84,29],[18,31],[49,31],[80,31],[19,41],[60,41],[13,43],[56,43],[1,53],[54,53],[46,61],[2,67],[6,67],[44,67],[69,67],[73,67],[50,73],[23,79],[47,79],[73,79],[28,89],[62,89],[73,89],

Aggiusta il factor base algebrico A={(r,p)} ordinandolo rispetto a p e con soli 20 elementi. Rimango in attesa

[[0,2],[6,7],[13,17],[11,23],[26,29],[18,31],[19,41],[13,43],[1,53],[46,61],[2,67],[6,67],[44,67],[50,73],[23,79],[47,79],[73,79],[28,89],[62,89],[73,89]];

FINE SETUP

== A ==
factor base algebrico A
[[0, 2], [6, 7], [13, 17], [11, 23], [26, 29], [18, 31], [19, 41], [13, 43], [1, 53], [46, 61], [2, 67], [6, 67], [44, 67], [50, 73], [23, 79], [47, 79], [73, 79], [28, 89], [62, 89], [73, 89]]

== Q ==
character base quadratico Q
List([[28, 97], [87, 101], [47, 103], [4, 107], [8, 107], [80, 107]])

Grado massimo del polinomio f: 3

Cerca coppie (a,b) smooth in R e A (primi ideali)

1: (-73 ,1 )             smooth in R e A
2: (-13 ,1 )             smooth in R e A
3: (-6 ,1 )              smooth in R e A
4: (-2 ,1 )              smooth in R e A
5: (-1 ,1 )              smooth in R e A
6: (1 ,1 )               smooth in R e A
7: (2 ,1 )               smooth in R e A
8: (3 ,1 )               smooth in R e A
9: (13 ,1 )              smooth in R e A
10: (15 ,1 )             smooth in R e A
11: (23 ,1 )             smooth in R e A
12: (61 ,1 )             smooth in R e A
13: (1 ,2 )              smooth in R e A
14: (3 ,2 )              smooth in R e A
15: (33 ,2 )             smooth in R e A
16: (2 ,3 )              smooth in R e A
17: (5 ,3 )              smooth in R e A
18: (19 ,4 )             smooth in R e A
19: (14 ,5 )             smooth in R e A
20: (37 ,5 )             smooth in R e A
21: (11 ,7 )             smooth in R e A
22: (15 ,7 )             smooth in R e A
23: (-7 ,9 )             smooth in R e A
24: (115 ,11 )           smooth in R e A
25: (119 ,11 )           smooth in R e A
26: (175 ,13 )           smooth in R e A
27: (5 ,17 )             smooth in R e A
28: (33 ,17 )            smooth in R e A
29: (-1 ,19 )            smooth in R e A
30: (35 ,19 )            smooth in R e A
31: (17 ,25 )            smooth in R e A
32: (49 ,26 )            smooth in R e A
33: (375 ,29 )           smooth in R e A
34: (9 ,32 )             smooth in R e A
35: (1 ,33 )             smooth in R e A
36: (5 ,41 )             smooth in R e A
37: (9 ,41 )             smooth in R e A

N.ro elementi trovati 37

Creazione matrice esponenti modulo 2

Qui vi dirà la matrice E.

Al momento ci arrestiamo qua e in un altro post vedremo la parte finale per la ricerca delle radici.

Sorgente Comb.gp

Lo trovate sui post di CFRAC e contiene la funzione CombList().

Sorgente SetupGNFS.txt

\\ SETUP iniziale
\\ n=45113
m=31;          \\ base espansione del polinomio
a1=-400;a2=400; \\ intervallo di a
b1=1;b2=41; \\ intervallo di b, che
            \\ non può essere a cavallo di zero
dimR=10;    \\ dimensione di R
lim_pR=30;  \\ limite max dei primi di R
lim_pA=90;  \\ limite max dei primi di A
lim_pQ=107;  \\ limite dei primi di Q. u=|Q|
                        \\ i primi di Q non devono essere in A
                        \\ l=|A|volte_AtoR=2; \\ la dimensione di A rispetto ad R

\\ polinomio scelto
f(x) = x^3 + 15 * x^2 + 29*x + 8
\\ Fine SETUP iniziale


Sorgente GNFS.gp

\\ Author: Rosario Turco


\\ INCLUDE
\\ di routine di utilità
\\ modificare solo il path per Comb.txt

\r c:/projpari/Comb.txt
\r c:/projpari/SUMMOD2.txt

\\ SETUP
\\ carico il Setup
\r c:/projpari/SetupGNFS.txt

\\ MEMO del file di setup
print("\n SimpleGNFS starta con la configurazione di SetupGNFS.txt - Verificare se adatta\n Per partire digita: simpleGNFS(n) Es. simpleGNFS(45113)\n ")

\\ SOFTWARE
\\ n numero da fattorizzare es.45113
simpleGNFS(n)=local(d=0);{

\\ VERIFICA se PRIMO
if(isprime(n)==1, print("\nNumero primo"); return;);
if(isprime(n)==0, print("\nNumero composto\n"););

print("Inizio SETUP");
print("\nPolinomio configurato in SETUP");

print("\nfactor base R di ",dimR," elementi");
R=listcreate(dimR);
i=1;
\\ ne prendo solo dimR con primi minore di lim_pR
forprime(p=1,10*dimR,if(i<=dimR & pprint("\n== R ==");
print("factor base R");
print(R);

if(length(R)

print("\nfactor base algebrico A. Esempio:");
forprime(p=1,lim_pA,for(r=0,lim_pA,if(f(r)%p==0,print1("[",r,",",p,"],");)));
A = listcreate(volte_AtoR*dimR);
print("\n\nAggiusta il factor base algebrico A={(r,p)} ordinandolo rispetto a p e con soli ",volte_AtoR*dimR," elementi. Rimango in attesa\n");
A=input();

print("\n== A ==");
print("factor base algebrico A");
print(A);

if(length(A)!=volte_AtoR*dimR,print("dim di A errata");return;);
if(lim_pQ<=lim_pA,print("Aumentare lim_pQ rispetto a lim_pA");return;);

contQ=0;
forprime(q=A[2*dimR][2]+1,lim_pQ,for(s=0,lim_pA,if(f(s)%q==0,contQ++)));
Q=listcreate(contQ);
forprime(q=A[2*dimR][2]+1,lim_pQ,for(s=0,lim_pA,if(f(s)%q==0,listput(Q,[s,q]);)));

print("\n== Q ==");
print("character base quadratico Q");
print(Q);

if(length(Q)>=length(R), print("dim Q >= dim R. Diminuire in SETUP lim_pQ o aumentare dimR"); return;);


d=maxm(n,m); \\ grado d del polinomio f
print("\nGrado massimo del polinomio f: ",d);

print("\nCerca coppie (a,b) smooth in R e A (primi ideali) \n");
cerca_coppie_ab(n,m,a1,a2,b1,b2,d);

print("\nCreazione matrice esponenti modulo 2\n");
 createE(m,d);
 \\ Cerca soluzione
print("\nCerca soluzione non banale\n");
   \\ SumEMod2(E,n);
 \\ Cerca gcd
 \\ Fine elaborazione
}

\\ ROUTINE

/*
**    maxm(n,m)
**    Ritorna m^k <= n
*/

maxm(n,m) = local (k);{
 k=0;
 while( m^k <= n, k++;);
 return(k-1);
}

/*
**    Cerca coppie (a,b) smooth in R e in A
**
*/
cerca_coppie_ab(n,m,a1,a2,b1,b2,d)=local(b=0,a=0,x=0,cont=0);{
  global(FinalSieving);
  FinalSieving=listcreate(100);

  for(b=b1,b2,
    for(a=a1,a2,
 x=a+b*m;
 if(is_inR(x)==1, \\ smooth in R
   \\se a+bT è divisibile allora
   \\ metti in lista (a,b)

if(is_div(abs(((b)^d) * f(-a/b)))==1,
   \\ smooth in R e in A
print(cont+1,": (",a," ,",b," ) \t\t smooth in R e A");
listput(FinalSieving,[a,b]);
cont++;
);

 );
 x=0;
);
  );
  print("\nN.ro elementi trovati ", cont);

<1 1="" arning="" di="" elementi="" k="" l="" length="" minore="" n.ro="" nbsp="" p="" print="" trovati="" u=""><1 1="" arning="" di="" elementi="" k="" l="" length="" minore="" n.ro="" nbsp="" p="" print="" trovati="" u=""><1 1="" arning="" br="" di="" elementi="" k="" l="" length="" minore="" n.ro="" nbsp="" print="" trovati="" u="">
<1 1="" arning="" di="" elementi="" k="" l="" length="" minore="" n.ro="" nbsp="" p="" print="" trovati="" u=""><1 1="" arning="" br="" di="" elementi="" k="" l="" length="" minore="" n.ro="" nbsp="" print="" trovati="" u="">  return(FinalSieving);
}

/*
**    Verifica se x è in R
*/
is_inR(x)=local(cont=0,i=1);{
  y=factor(abs(x));
  for(k=1,matsize(y)[1],
    while(i<=length(R),
 if( y[k,1]==R[i],  
     cont++
 );
 i++;
    );
    i=1;
  );
  if(cont==matsize(y)[1], return(1););
  return(0);
}

/*
** is_div
** deve verificare se tra tutte le combinazioni Clist
** di prodotti di elementi in A ce n'è uno che vale
** quanto value
** Se esiste is_div ritorna 1 altrimenti ritorna 0
*/
is_div(value)=local(pr=1);{
  global(Clist);
  CombList(length(A));
  for(k=1,length(Clist),
    for(i=1,length(Clist[k]),
       pr=pr*A[Clist[k][i]][2];
);
if(pr==value,
 return(1);
);
pr=1;
  );
  return(0);
}

/*
**  Crea la matrice modulo 2 degli esponenti
*/
createE(m,d)=local();{
  global(FinalSieving,A,R,Q);
  col=1+length(R)+length(A)+length(Q); E=matrix(length(FinalSieving),col);

  for(i=1,length(FinalSieving),
    \\ segno 1 o 0
x= FinalSieving[i][1]+m*FinalSieving[i][2];
if(x>=0,E[i,1]=0;);
    if(x<0 i="" p=""> 
\\ parte di R
    w=factor(x);
for(j=1, matsize(w)[1],
 pos=SearchInR(w[j,1]);
 if(pos!=0,
   E[i,pos+1]=w[j,2]%2;
 );
 pos=0;
);
\\ Parte di A
a=FinalSieving[i][1];
b=FinalSieving[i][2];
x=abs(((b)^d) * f(-a/b));
w=factor(x);
for(j=1, matsize(w)[1],
 pos=SearchInA(2,w[j,1],volte_AtoR*dimR);
 if(pos!=0,
    E[i,1+length(R)+pos]=w[j,2]%2;
 );
 pos=0;
);
    \\ parte di Q
for(k=1,length(Q),
   \\ simbolo di Legendre
   if(kronecker(a+b*Q[k][1],Q[k][2])==1,
    E[i,1+length(R)+length(A)+k]=0;
);
\\ simbolo di Legendre
        if(kronecker(a+b*Q[k][1],Q[k][2])!=1,
     E[i,1+length(R)+length(A)+k]=1;
);
);
  );
 return(E);
}

SearchInR(ele)=local();{
  for(i=1,length(R),
if(ele==R[i],return(i););
  );
  return(0);
}

SearchInA(start,ele,dim)=local(ret=0);{
  \\ ritorna ultima posizione con quel valore
  for(i=start,dim,
if(ele==A[i][2],ret=i;);
  );
  return(ret);
}

Sorgente SUMMOD2.gp
Questo sorgente lo vedremo sviluppato meglio in un prossimo post. Ed è la chiave per passare dalla matrice alla soluzione.

<1 1="" arning="" di="" elementi="" k="" l="" length="" minore="" n.ro="" nbsp="" p="" print="" trovati="" u=""><1 1="" arning="" di="" elementi="" k="" l="" length="" minore="" n.ro="" nbsp="" p="" print="" trovati="" u=""><1 1="" arning="" br="" di="" elementi="" k="" l="" length="" minore="" n.ro="" nbsp="" print="" trovati="" u=""><0 i="" p="">
<1 1="" arning="" di="" elementi="" k="" l="" length="" minore="" n.ro="" nbsp="" p="" print="" trovati="" u=""><1 1="" arning="" di="" elementi="" k="" l="" length="" minore="" n.ro="" nbsp="" p="" print="" trovati="" u=""><1 1="" arning="" br="" di="" elementi="" k="" l="" length="" minore="" n.ro="" nbsp="" print="" trovati="" u=""><0 i="" p="">Alla prox





domenica 22 gennaio 2017

GNFS - Esempio - parte 3



Siamo arrivati all'esempio. Giuro! 

Discuteremo pero l'esempio in modo da consolidare i concetti visti nei post precedenti, parte 1 e parte 2, specie se vorrete automatizzare il tutto con del software.

Esempio n=45113.

Faremo prima di tutto un test per stabilire se n risulta composto. Ad esempio in PARI GP:
gp > n=45113;
gp > isprime(n)
%2 = 0
quindi composto. Proviamo a fattorizzare col GNFS.

I primi 3 step sono di setup; il setup lo cerchiamo di dimensionare in modo tale da non avere un grosso tempo di elaborazione in fase di algebra lineare. Però facendo in questo modo non vuol dire che riusciremo a fattorizzare subito e dovremmo, poi, ricominciare daccapo il setup, allargando le dimensioni.

Step 1: selezionare m e f. 
Supponiamo arbitrariamente m=31.Il modo migliore per ottenere una f ed un f(m) = 0 mod n comporta di espandere n in base m. 

45113 = 31^3 + 15* 31^2 + 29*31 + 8

Quindi il polinomio f scelto risulta:
f(x) = x^3 + 15 * x^2 + 29*x + 8

da cui

f(31) = 45113 
f(31) = 0 mod 45113

Step 2: trovare fattore base razionale R e fattore base algebrico A.

Indicheremo nel seguito col simbolo |S| la dimensione dell'insieme S generico.

Una considerazione. La dimensione dei fattor base viene scelta arbitrariamente, ma più sono grandi, maggiore sara il tempo di elaborazione che si otterrà nello step di algebra lineare. Per cui inizialmente non conviene esagerare e si aumenta in seguito solo se nei tentativi fatti, di impostazione del GNFS,non si è riusciti a fattorizzare.

Per il fattore base razionale R inizialmente consideriamo {pi} con pi<30 div="">
gp > R=listcreate(10);
gp > i=1;
forprime(p=1,100,if(i<=10 & p<30,listput(R,p);i++));
R
%12 = List([2, 3, 5, 7, 11, 13, 17, 19, 23, 29])

Quindi:
R = [2,3,5,7,11,13,17,19,23,29], |R|=k=10

Per il fattore base algebrico A, ricordando il Teorema 3, dobbiamo scegliere delle coppie (ri,p) dove p e primo e f(ri) = 0 mod p.

Scegliamo, anche qui arbitrariamente, <90 cercando="" che="" cio="" di="" fare="" i="" in="" modo="">la dimensione di A che sia circa 2-3 volte maggiore di quella di R cioe l>k. )iniziamo con 2).


Cerchiamo ora un insieme {ri} tale che, per ogni i, risulti f(ri)=0 mod p. Per ogni ri trovato aggiungiamo, poi, la coppia (ri,p) ad A. Per la considerazione fatta prima, non necessariamente l'insieme di ri deve essere tutto quello che si può ottenere. 

Dovremmo pero ripetere il procedimento per tutti i p minori di 90, quindi possiamo usare PARI GP per creare la lista A:
(17:33) gp > f(x) = x^3 + 1*x^2 + 29*x + 8;
(17:33) gp > (17:50) gp > forprime(p=1,90,for(r=0,90,if(f(r)%p==0,print1("[",r,",",p,"],"))));
[0,2],[2,2],[4,2],[6,2],[8,2],[10,2],[12,2],[14,2],[16,2],[18,2],[20,2],[22,2],[24,2],[26,2],[28,2],[30,2],[32,2],[34,2],[36,2],[38,2],[40,2],[42,2],[44,2],[46,2],[48,2],[50,2],[52,2],[54,2],[56,2],[58,2],[60,2],[62,2],[64,2],[66,2],[68,2],[70,2],[72,2],[74,2],[76,2],[78,2],[80,2],[82,2],[84,2],[86,2],[88,2],[90,2],[1,3],[2,3],[4,3],[5,3],[7,3],[8,3],[10,3],[11,3],[13,3],[14,3],[16,3],[17,3],[19,3],[20,3],[22,3],[23,3],[25,3],[26,3],[28,3],[29,3],[31,3],[32,3],[34,3],[35,3],[37,3],[38,3],[40,3],[41,3],[43,3],[44,3],[46,3],[47,3],[49,3],[50,3],[52,3],[53,3],[55,3],[56,3],[58,3],[59,3],[61,3],[62,3],[64,3],[65,3],[67,3],[68,3],[70,3],[71,3],[73,3],[74,3],[76,3],[77,3],[79,3],[80,3],[82,3],[83,3],[85,3],[86,3],[88,3],[89,3],[6,7],[13,7],[20,7],[27,7],[34,7],[41,7],[48,7],[55,7],[62,7],[69,7],[76,7],[83,7],[90,7],[1,13],[2,13],[9,13],[14,13],[15,13],[22,13],[27,13],[28,13],[35,13],[40,13],[41,13],[48,13],[53,13],[54,13],[61,13],[66,13],[67,13],[74,13],[79,13],[80,13],[87,13],[4,17],[8,17],[21,17],[25,17],[38,17],[42,17],[55,17],[59,17],[72,17],[76,17],[89,17],[6,31],[37,31],[68,31],[25,37],[62,37],[14,43],[57,43],[34,47],[81,47],[44,53],[24,59],[83,59],[31,61],[7,67],[16,67],[43,67],[74,67],[83,67],[58,73],[74,79],[9,83],[35,89],

Togliamo l'ultima virgola e definiamo la lista. Scegliamo una size |A|=l=20 per cui da tutte queste coppie (r,p) ne scegliamo solo 20, per cui A diventa:

A = listcreate(20); A=[[0,2],[6,7],[13,17],[11,23],[26,29],[18,31],[19,41],[13,43],[1,53],[46,61],[2,67],[6,67],[44,67],[50,73],[23,79],[47,79],[73,79],[28,89],[62,89],[73,89]]

Ora dobbiamo scegliere un fattore base quadratico Q.
Scegliamo numeri primi qi che non sono presenti in A e per ogni qi scelto cerchiamo degli si tali che f(si)=0 mod qi. Poi, per ogni si trovato, aggiungiamo a Q la coppia (si,qi).

Prendiamo i numeri primi qi=97,101,103 e 107 non presenti in A. In PARI GP:
gp > forprime(q=97,107,for(s=0,90,if(f(s)%q==0,print1("[",s,",",q,"],"))));
[28,97],[87,101],[47,103],[4,107],[8,107],[80,107],

Togliamo anche qui l'ultima virgola e definiamo la lista.

Quindi 
Q=listcreate(6);
Q=[[28,97],[87,101],[47,103],[4,107],[8,107],[80,107]]
|Q|=u=6

Le coppie (a,b) da cercare saranno, quindi, in totale 1+k+l+u=37 con a+bm smooth in R e a+bT smooth in A.

Step 3: creare i due vettori o liste di sieving.

Decidiamo, arbitrariamente anche qui, che b sia nell'intervallo (1,41) mentre a sia nell'intervallo (-400,400). In realtà dopo un paio di setup si trova il valore dell'intervallo di a e di b da usare come anche le dimensioni di R, A e Q.

Dovremo ciclare sull'intervallo di possibili valori di b, e per ogni b, costruire due liste, ognuna con  (2*400+1) elementi. 

Con le due liste dobbiamo cercare, poi, quei valori di a, per i quali a+bm e a+bT sono smooth su R ed A contemporaneamente. Se li troviamo, salviamo la coppia (a,b) in una lista finale che potremmo chiamare ListSmoothFinalSieving. 

Ricordiamo che un numero intero I in Z e detto smooth (liscio) sopra un fattore base B se B contiene tutti i fattori primi di I.

Ad esempio per la coppia (119,11) si puo vedere che a+bm=119+11*31=22*5*23 e questi fattori sono tutti compresi in R, quindi la coppia (119,11) e smooth su R.

Le seguenti coppie (r,p) in A, che abbiamo trovato precedentemente, dividono a+bT = 119+11*T: {(19,41),(44,67),(62,89)}. 

Difatti dal Teorema 5, questa è una completa fattorizzazione di a+bT se e solo se 

gp > 41*67*89 == abs((-11)^3 * f(-119/11))

Uscendo 1 significa che è vero.

Ricaviamoci ora ListSmoothFinalSieving. Per fare questo ci conviene scrivere del software, cerca_coppie_ab, che farò vedere in un altro post ma che segue tutta la teoria vista finora.

Supponiamo al momento che la lista ListSmoothFinalSieving contenga i seguenti valori:

%217 = List([[-73, 1], [-13, 1], [-6, 1], [-2, 1], [-1, 1], [1, 1], [2, 1], [3,1], [13, 1], [15, 1], [23, 1], [61, 1], [1, 2], [3, 2], [33, 2], [2, 3], [5, 3], [19, 4], [14, 5], [37, 5], [11, 7], [15, 7], [-7, 9], [115, 11], [119, 11], [175, 13], [5, 17], [33, 17], [-1, 19], [35, 19], [17, 25], [49, 26], [375, 29], [9, 32], [1, 33], [5, 41], [9, 41]])

Step 4: creare la matrice X (algebra lineare)

La matrice X prodotta con un insieme di vettori V. Ad esempio il vettore V corrispondente alla coppia (119,11) corrisponde a:
- primo elemento 0 (segno di a+bm)
- esponenti modulo 2 dei fattori di a+bm: 0010000010   (esponenti riferiti ai primi di R)
- esponenti modulo 2 dei fattori a+bT: 00000010000010000010 (esponenti riferiti ai primi di A)
- esponenti modulo 2 dei fattori di Q: 110000  (esponenti riferiti ai primi di Q)

Una possibile soluzione dell'equazione

[A1,...,Ay]^T = [0,0,0,1,0,1,0,0,1,1,1,0,0,1,1,0,1,1,1,0,0,0,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0]^T

Questo implica che per le seguenti coppie (a,b)
V = {(−2,1),(1,1),(13,1),(15,1),(23,1),(3,2),(33,2),(5,3), (19,4),(14,5),(15,7),(119,11),(175,13),(−1,19),(49,26)}

\[a) \prod_{{(a,b)} \in V} (a+bm)\] è un quadrato perfetto in Z e che 
\[ b) \prod_{{(a,b)} \in V} (a+bT)\] è un perfetto quadrato in Z[T]. 

Qui 
a) = 45999712751795195582606376960000 
b) = 58251363820606365·T^2 +149816899035790332·T
 +75158930297695972 
 A questo punto occorre computare le radici. 
 E' facile notare che:
 a) = 2553045317222400^2
 b) = (108141021·T^2 +235698019·T +62585630)^2
Usando il mapping del Teorema 2, 
φ(108141021·T^2 +235698019·T +62585630)= 111292745400 

Quindi 1112927454002≡2553045317222400^2 (mod n)

Step 5: cerchiamo il gcd
Da qui
gcd(45113,111292745400+2553045317222400) = 197 gcd(45113,111292745400−2553045317222400) = 229 

Per cui n = 45113 = 197·229. Il GNFS quindi è riuscito a fattorizzare.

Se avrò tempo vi farò vedere in PARI GP il mio SimpleGNFS!

Alla prox

sabato 14 gennaio 2017

GNFS, teoria - parte 2



Ancora un po' di teoria e chiarimenti, in questa parte 2, per affinare ancora qualcosa prima degli esempi. Mi odierete lo so, ma un poco di pazienza e ci siamo quasi all'esempio. Anche perchè non solo ho dovuto pensare agli approfondimenti ma a mettervi le formule in LaTex.

Precedentemente abbiamo visto che dobbiamo cercare un insieme finito di numeri U ={(a,b)} tali che a+bm sia smooth su un fattore base razionale R e a+bT sia smooth su un fattore base algebrico A.

Sia R costituito da k elementi, A con l elementi. Scegliamo un arbitrario carattere base quadratico Q con u elementi.

R e A li useremo per cercare un quadrato in Z e Z[T], mentre Q sarà usato per verificare che il risultato sia un quadrato.

Ogni coppia (a,b) in U può essere rappresentata come un vettore riga con 1+k+l+u ingressi (o entry).

Il primo ingresso del vettore riga deve essere 0 se a+bm è positivo, oppure vale 1 se a+bm e negativo.

I successivi k ingressi sono il vettore degli esponenti modulo 2.

I successivi l elementi indicano se un particolare elemento di A divide a + bT. Ovvero qui appaiono gli esponenti  della fattorizzazione modulo 2.

Gli ultimi u ingressi sono usati in connessione con Q. Ogni entry è settato a 0 se per la relativa coppia (s,q) e:
\[\big(\frac{a+bs}{q} \big) = 1\] Altrimenti viene settato a 1. Quindi sia:
- R = {t1,t2,...,tk}
- A = {(r1,p1),(r2,p2),...(rl,pl)}
- Q = {(s1,q1),(s2,q2),...(su,qu)}.

Per una data coppia (a,b), a+bm ha una fattorizzazione
\[t_{1}^{e_{1}} t_{2}^{e_{2}} ... t_{k}^{e_{k}}.\] a+bT ha una fattorizzazione che puo essere rappresentata da:
\[(r_{1},p_{1})^{f_{1}}(r_{2},p_{2})^{f_{2}} ...(r_{l},p_{l})^{f_{l}}.\]A questo punto la coppia (a,b) si puo rappresentare con un vettore riga del tipo:
- primo elemento: 0 se a+bm>=0 oppure 1;
- k elementi: e1(mod2),e2(mod2),...,ek(mod2)
- l elementi: f1(mod2),f2(mod2),...,fl(mod2)
- u elementi: per ogni i=1...u, e 0 se \[\big(\frac{a+bs_{i}}{q_{i}} \big) = 1\] altrimenti vale 1.

Dobbiamo giustificare ora un attimo come siamo arrivati a questo vettore.

Supponiamo ora che un V incluso in U sia trovato tale che siano quadrati perfetti in Z e Z[T] entrambi :\[ a) \prod_{{(a,b)} \in V} (a+bm)\] \[ b) \prod_{{(a,b)} \in V} (a+bT)\] allora ne conseguono diverse cose:

1) il produttorio a) deve essere positivo. Indichiamo con cj il valore 0 se a+bm>=0 oppure 1, cioè il primo ingresso del vettore. Allora a) è positivo se e solo se \[\sum c_{j} = 0(mod2) \] Questo assicura che il numero di numeri negativi nel prodotto è pari. Inoltre -1 elevato ad un pari da 1.

2) ogni esponente della fattorizzazione di a) deve essere pari. Questo perchè per ogni i, aj+bjm è smooth su R e lo sarà anche il prodotto. Inoltre se ej è l'esponente che appare per qualche t appartenete ad R nella fattorizzazione di a) allora tale esponente apparirà anche nel prodotto per cui la a) è quadrata per ogni t appartenente ad R. In altri termini e:\[\sum_{(a_{j},b_{j}) \in V} e_{j} = 0  \Leftrightarrow  \sum_{(a_{j},b_{j}) \in V} (e_{j} mod 2)  \equiv  0 (mod2)\] Ogni ej mod2 e un ingresso del vettore.

3) ogni esponente che si ha nella fattorizzazione della b) deve essere pari. Come prima, per ogni i, aj+bjT è smooth su A e lo sara anche il prodotto. Poichè aj+bjT hanno la rappresentazione
(r1,p1)^ej1(r2,p2)^ej2 ... (rl,pl)^ejl allora la b) ha per l elementi una rappresentazione \[(r_{1},p_{1})^{\sum e_{j_{1}}}(r_{l},p_{l})^{\sum e_{j_{l}}} \] Ogni esponente sarà qui pari per gli stessi ragionamenti del punto precedente per ogni i. Tale cosa rappresenta l'entry del vettore

4) Per ogni (s,q) in Q, \[\prod_{{(a_{j},b_{j})} \in V} \big(\frac{a_{j}+b_{j}s}{q} \big) = 1\] A causa del fatto che per ogni j il simbolo di Legendre deve appartenere all'intervallo (-1,1), allora per fare in modo che il precedente produttorio sia uguale a 1 occorre che il numero di j, tali che il simbolo di Legendre è uguale a -1, sia pari. Per una data coppia (s,q) la rappresentazione del vettore di (aj,bj) ha una entry del tipo vale 0 se il simbolo di Legendre (aj+bjs/q)=1 oppure vale 1 altrimenti. Se la somma di queste entry è pari allora il numero dei j elementi per i quali il simbolo di Legendre è uguale a 1 è pari. Quindi:\[\prod_{{(a_{j},b_{j})} \in V}  \big(\frac{a_{j}+b_{j}s}{q} \big) = 1 \Leftrightarrow  \sum_{(a_{j},b_{j}) \in V} {\begin{cases}0,\prod_{{(a_{j},b_{j})} \in V} \big(\frac{a_{j}+b_{j}s}{q} \big) = 1 \\1,else\end{cases}}\]

Poichè devono essere vere contemporaneamente sia 1),2),3),4) allora affinchè \[\prod_{{(a_{j},b_{j})} \in V}  (a_{j}+b_{j}m)\] è un quadrato perfetto in Z e \[\prod_{{(a_{j},b_{j})} \in V}  (a_{j}+b_{j}T)\] è un quadrato perfetto in Z[T] se e solo se la somma modulo 2 dei vettori che sono rappresentazione di ogni (aj,bj) danno come risultato un vettore con tutti 0.

Per dire quest'ultimo fatto in modo più rigoroso, sia U un insieme di numeri smooth (a,b) avente y elementi. Sia X una matrice di y*(1+k+l+u) elementi dove ogni riga è un vettore della rappresentazione di (a,b) in U. Ora cercare un V incluso in U per ottenere un quadrato perfetto è un passo di algebra lineare cioè si deve cercare un vettore colonna A tale che: \[X^{T}\begin{bmatrix}A_{1} \\A_{2}\\...\\A_{y} \end{bmatrix} \equiv 0 (mod 2)\]

Se y>1+k+l+u allora la congruenza ha sicuramente una soluzione non banale.

L'algoritmo GNFS
Definiamo qui l'algoritmo del GNFS.

Il GNFS per fattorizzare un numero composto n generico fa i seguenti passi:

a) Si sceglie un m appartenete a Z e si cerca un corrispondente f tale che f(m) = 0 (mod n) sfruttando il metodo dell'espansione base m del numero.

b) Si definisce un fattore base razionale R tale che R ha finiti e molti elementi e per ogni x in R, x è prime. Sia k il numero di elementi in R.

c) Si definisce un fattore base algebrico A tale che A è un numero finito di molti elementi e per ogni (r,p) in A, p è primo e r soddisfa f(r)=0 (mod p). Sia l il numero di elementi in A.

d) Si definisce un carattere base quadratico Q con finiti molti elementi tali che per ogni (s,q) in Q, q è primo e f(s) = 0 (mod q). Assicurarsi che per ogni (s,q) in Q, (s,q) non appartenga ad A. Sia u il numero di elementi in Q.

e) Per un fissato b in Z, costruiamo un vattore di sieving come visto nella parte 1. Osserviamo gli elementi del vettore di sieving che sono smooth e per quegli elementi registriamo quali qi in R e quali (ri,pi) in A sono divisori. Ripetiamo questo processo per vari b finchè abbiamo più di 1+k+l+u coppie(a,b) trovate tale che a + bm sia smooth in Z e a + bT e smooth in Z[T]. Sia y il numero di coppie smooth (a,b) trovate.

f) Popoliamo una matrice X di y * (1+k + l + m) righe per colonne X come descritto prima.

g) Risolviamo l'equazione
\[X^{T}\begin{bmatrix}A_{1} \\A_{2}\\...\\A_{y} \end{bmatrix} \equiv 0 (mod 2)\] per {A1,...,Ay}.
Sia V incluso in U, è un sottoinsieme definito da ogni (aj,bj) in U,(aj,bj) appartenente a V se Aj = 1. Allora \[\prod_{{(a_{j},b_{j})} \in V}  (a_{j}+b_{j}m)\] è un quadrato perfetto in Z e \[\prod_{{(a_{j},b_{j})} \in V}  (a_{j}+b_{j}T)\] è un quadrato perfetto in Z[T].

h) Con il mapping visto nella parte 1
\[\phi(\prod_{a_{j},b_{j} \in V} (a+bT)) \equiv \prod_{a_{j},b_{j} \in V} (a_{j}+b_{j}m)(mod2)\]

Si usa ciò per arrivare a fattorizzare  n con la differenza di quadrati. Se nessuna fattorizzazione e trovata si ritorna a b) e si ripete il processo.

Nella parte 3 faremo finalmente l'esempio.

Alla prox




lunedì 9 gennaio 2017

GNFS, teoria - Parte 1

Il GNFS (General Number Field Sieve) si basa su una teoria evoluta e abbastanza consolidata. Proverò a ridurre il tutto alle cose essenziali.

Il GNFS discende da una serie di metodi che sono passati attraverso varie innovazioni principali: metodo di Fermat, SQUFOF, CFRAC, Quadratic Sieve, Number Field Sieve. In questo percorso, dal Settecento agli anni '80 del Novecento, poi sono stati introdotte delle varianti, più che altro dei tecnicismi e ottimizzazioni prestazionali su parti del metodo: ad esempio si è passati dal metodo di eliminazione di Gauss (di algebra lineare) a quello a blocchi di Lanczos per trattare matrici sparse enormi,oppure si è passati al parallelismo, interno ai task nella stessa fase. Non è però possibile parallelizzare le fasi tra loro, ovvero il risultato di una fase è propedeutico alla successiva.

Sulla fattorizzazione, quindi, si è fermi ai risultati del Novecento e nessuna nuova e importante innovazione è finora nata.

Il GNFS allo stato attuale è il miglior algoritmo in termini di velocità e complessità, finora noto.

Qui esamineremo il metodo base concettuale. Ci aiuteranno anche le piccole somiglianze col CFRAC, esaminato nei tre post precedenti.Vi darò anche dei riferimenti per l'approfondimento.

Col GNFS, per fattorizzare un intero n, abbiamo bisogno di alcune cose.

Innanzitutto ci servono:
- un polinomio f: R -> R, a coefficienti interi e irriducibile (Vedi [1] per definizione di polinomio irriducibile).
- un parametro m tale che f(m)= 0 mod n
- un anello Z[T] su cui sono definiti i polinomi

Approfondiamo un pò le cose di cui sopra.

Solitamente m si sceglie facendo l'espansione del numero n. Se \[n = a_{d}*x^{d} + a_{d-1}*x^{d-1} + ... + a_{0}\]
allora si sceglie m in modo che f(m) = n  e quindi f(m) = 0 mod n. Una buona tecnica (vedi [3]) per avere polinomi di qualità è di scegliere polinomi con coefficienti interi piccoli.

Sia \[T \in \mathbb{C}\] una radice del polinomio f di cui sopra, cioè tale che \[f(T) = 0\] Lo spazio dell'anello sia ora definito come segue:

\[Z[T]=set x: x=a_{d}*x^{d} + a_{d-1}*x^{d-1} + ... + a_{0} for {a_{j} \subset  \mathbb{Z}}\] Teorema 1
Lo spazio Z[T] con l'operazione di moltiplicazione è un anello.

L'anello dei polinomi appena definito, a causa della moltiplicazione, ha delle particolari proprietà. Ora vediamo quali.

Siano A,B in Z[T] e a(x),b(x) due suoi polinomi, tali che a(T)=A e b(T)=B. Inoltre c(T)=C.

Per l'algoritmo della divisione, sia che
a(x)b(x)=e(x)f(x)+c(x), allora:

AB = a(x)b(x)|x=T
   = e(x)f(x)|x=T + c(x)|x=T
   = e(T)f(T)+ C
   = C

Questo poichè f(T)=0. Per costruzione, inoltre, il grado del polinomio c è minore di f (abbiamo fatto la divisione precedentemente).
 
Il niocciolo dell'algoritmo GNFS
Siano assegnati due quadrati perfetti: \[{\beta^{2} \in Z[T], y^{2} \in \mathbb{Z}}\]
Il tutto gira attorno al seguente

Teorema 2:
Sia f(x) un polinomio a coefficienti interi, sia \[T \in \mathbb{C}\] una radice di f, m in Z/nZ tale che f(m) = 0 mod n, allora esiste un solo mapping \[\phi: Z[T] -> Z/nZ\] che soddisfa le proprietà seguenti:

\[a) \phi(ab) = \phi(a)\phi(b), con a,b \in \mathbb{Z[T]}\]
\[b) \phi(a+b) = \phi(a) + \phi(b)\]
\[c) \phi(1) = 1 mod n\]
\[d) \phi(T) = m mod n\]
Che ci facciamo? Si può ottenere una differenza di quadrati (con i polinomi), un pò come col metodo di Fermat, quando x^2-y^2=(x+a)(x-a) e si arrivava a fattorizzare. Solo che qui occorre qualche passo in più con qualche sofisticazione in più, ma col vantaggio di fattorizzare numeri con moltissime cifre.

Supponiamo che esista un insieme U di coppie di interi (a,b) e beta in Z[T] tali che:

\[\prod_{{(a,b)} \in U} (a+bT) =  \beta^{2}  \]
\[\prod_{{(a,b)} \in U} (a+bm) =  y^{2}  \]

Se
\[x=\phi(\beta)\]
allora
\[x^{2}=\phi(\beta)\phi(\beta) = \phi(\beta^{2})\]

sostituendo  si ottiene che

\[= \phi(\prod_{{(a,b)} \in U}  (a+bT)) \]

per le proprietà del Teorema precedente si ha che:
\[= \prod_{{(a,b)} \in U} \phi(a+bT) \]
\[= \prod_{{(a,b)} \in U} (a+bm)\]
\[= y^{2}\]

Così, come facevamo nel CFRAC, riusciamo a trovare una relazione x^2 = y^2 mod n che ha 2/3 probabilità di ammettere una soluzione non banale e quindi fattorizzare!

Il punto di arrivo quindi è trovare quadrati perfetti in Z(T) e in Z e per farlo occorre seguire una strategia che viene inglobata nell'algoritmo.

Tale strategia, lunga, si basa su vari elementi necessari:
- fattore base razionale R e numeri smooth su R
- fattore base algebrico A e numeri smooth su A
- carattere base quadratico Q e carattere quadratico su Q
- algebra lineare, per lavorare sulle matrici
- gcd

Alcune cose ci sono note, altre le definiamo strada facendo, per non appesantire il tutto fin dall'inizio.

Se l'obiettivo è trovare dei quadrati perfetti in Z,
questo lo sappiamo fare da CFRAC. Supponiamo di conoscere un insieme di numeri i cui fattori primi non superano 19 (stiamo alludendo a numeri 19-smooth):
{455,39270,770,429,1616615,3990,106590,187,19019}

Una verifica immediata è vedere i fattori primi di ognuno con PARI GP
gp > L=[455,39270,770,429,1616615,3990,106590,187,19019]
gp > for(i=1,length(L), print(L[i],": \t", factor(L[i]));)
455:     [5, 1; 7, 1; 13, 1]
39270:   [2, 1; 3, 1; 5, 1; 7, 1; 11, 1; 17, 1]
770:     [2, 1; 5, 1; 7, 1; 11, 1]
429:     [3, 1; 11, 1; 13, 1]
1616615: [5, 1; 7, 1; 11, 1; 13, 1; 17, 1; 19, 1]
3990:   [2, 1; 3, 1; 5, 1; 7, 1; 19, 1]
106590: [2, 1; 3, 1; 5, 1; 11, 1; 17, 1; 19, 1]
187:     [11, 1; 17, 1]
19019:   [7, 1; 11, 1; 13, 1; 19, 1]

quindi un factor base è:
[2,3,5,7,11,13,17,19]

Inoltre tutti gli esponenti dei fattori sono uguali a 1.

Possiamo, rispetto al factor base, già costruire una prima matrice E degli esponenti modulo 2 :

2 3 5 7 11 13 17 19
-----------------------
0 0 1 1  0  1  0  0
1 1 1 1  1  0  1  0
1 0 1 1  1  0  1  0
0 1 0 0  1  1  0  0
0 0 1 1  1  1  1  1
1 1 1 1  0  0  0  1
1 1 1 0  1  0  1  1
0 0 0 0  1  0  1  0
0 0 0 1  1  1  0  1

Le combinazioni linearmente dipendenti sono le somme di vettori (o le righe della matrice E) per colonne che danno un vettore totale con tutti 0.

Se si gioca sulle combinazioni lineari dipendenti, è possibile trovare un fattore non banale come abbiamo visto in CFRAC.

Ad esempio se facciamo riferimento al prodotto
455·39270·1616615·3990, esso proviene dalla somma modulo 2 per colonna delle righe 1,2,5,6 della matrice E precedente, che dà come risultato un vettore totale di 0:

0 0 1 1  0  1  0  0
1 1 1 1  1  0  1  0
0 0 1 1  1  1  1  1
1 1 1 1  0  0  0  1
--------------------
0 0 0 0  0  0  0  0

In particolare, un semplice calcolo ci porta a dire che siamo di fronte a un quadrato perfetto:
 455·39270·1616615·3990=(339489150)^2.

Tutto questo ci ha insegnato che ci serve un fattore base su Z. Ma su Z[T] ci serve, invece, un fattore base razionale ed un fattore base algebrico.

La strategia richiede vari passi.

Step 1) Definizione di Smoothness su Z[T] e su Z
Un fattore base razionale R è un insieme finito di numeri primi. Un numero intero I \in Z è detto smooth (liscio) sopra un fattore base razionale R se R contiene tutti i fattori primi di I.

Se ritorniamo all'esempio di prima, il fattore base R contiene tutti i fattori primi dell'insieme dei numeri L= {455,39270,770,429,1616615,3990,106590,187,19019}

Un fattore base algebrico è un insieme finito {a+bT} incluso in Z[T] dove per a,b in Z,
ogni a+bT soddisfa la condizione che
\[\forall (a,b),\not\exists c,d  \in Z[T]\] tale che
c·d = a + bT. Questa condizione causa che a+bT sia chiamato "primo ideale".

Un elemento l in Z[T] è smooth sopra un fattore base algebrico A se \[\exists W \subset A: \exists  \prod_{{(c,d)} \in W} (c+dT) = l\] Questa definizione di fattore base algebrico però è poco manegevole, per cui ci riferiremo, come di solito si fa, ad un concetto equivalente dato dal Teorema successivo.

Teorema 3
Sia f(x) un polinomio a coefficienti interi, con T in C radice di f. Allora un insieme di coppie {(r,p)}, dove p è un numero primo ed r in Z/nZ con f(r)= 0 mod n, è in corrispondenza biunivoca con a+bT in Z[T] che soddisfa il criterio di essere un fattore base algebrico.

Cioè il Teorema 3 ci dice che al posto del fattore base algebrico a+bT possiamo usare un insieme finito di coppie di interi (r,p), che è equivalente.Tuttavia non tutti gli elementi di Z[T] possono essere rappresentati come coppie (r,p), però quello che può essere rappresentato torna utile al GNFS.

Step 2) Tecnica di setaccio (Sieving)

Sia R un arbitrario fattore base razionale rappresentato da un insieme di primi {qi} e sia A un arbitrario fattore base algebrico \in Z[T] rappresentato da un insieme di coppie {(ri,pi)} come desscritto dal Teorema 3.

Abbiamo bisogno di almeno altri 3 teoremi, per proseguire.

Teorema 4
Per c + dT, in un fattore base algebrico rappresentato da (r,p), c + dT divide a + bT in Z[T] se e solo se a = -br (mod p). 

Teorema 5 
Un insieme finito U di coppie (r,p) in Z[T] rappresentano una fattorizzazione completa di a+bT se e solo se 
\[\prod_{{(ri,pi)} \in U} p_{i} = (-b)^{d} * f(-a/b)\] dove d è il grado di f.

Teorema 6
Un numero primo q dividerà a+bm se e solo se a =-bm (mod q).

Con i teoremi 4,5,6 ora siamo pronti a trovare elementi smooth in Z[T] e Z con i seguenti passi di setaccio:
a) fissiamo b in Z e sia N un intero positivo.
b) creiamo due vettori colonna di setaccio, uno per le variazioni dei valori a + bT ed un altro per le variazioni dei valori a + bm, facendo variare a da -N a N: 

Primo vettore di setaccio
-N + bT
(-N+1) + bT
(-N+2) + bT
.................
N + bT

Secondo vettore di setaccio
-N + bm
(-N+1) + bm
(-N+2) + bm
.................
N + bm

c) Per ogni qi in R, sappiamo che qi dividerà a + bm se e solo se a =-bm (mod qi). Cercheremo, quindi, valori a per i quali a =-bm+kqi per qualche k \in Z e per ogni valore di a ci annoteremo il valore di a+bm nel secondo vettore di setaccio. Ripeteremo questo processo per ogni qi di R. Appena finito, prenderemo nota di tutti gli a+bm del vettore che sono fattori base con questo metodo. Questi a+bm sono smooth in R.

d) Usiamo un procedimento analogo per a+bT. Sappiamo che un (ri,pi) in A divide a + bT se e solo se a =-bri (mod pi). Quindi cerchiamo valori a tali che a =-bri+kpi per qualche k in Z. 
Per ogni a trovato, prendiamo nota della coppia (ri,pi) fattore base con a+bT nel primo vettore di setaccio. 

Appena finito, per tutti gli a+bT nel setaccio vi sarà una lista di fattori (ri,pi). 

Se \[\prod p_{i} = (-b)^{d} * f(-a/b)\] allora questa è una lista di elementi del fattore base e a+bT è smooth sul fattore base algebrico A.

e)Per ogni entry dobbiamo confrontare i due vettori. Se entrambi a+bm e a+bT sono smooth allora abbiamo trovato una relazione-coppia (a,b) da conservarci per il procedimento successivo.

Tutto il procedimento può essere ripetuto variando b e ottenendo le coppie (a,b).

Step 3) Verificare che gli elementi in Z[T] e Z sono quadrati. 

Questa è la parte un pò più lunga da digerire ma la più importante. Testare un s in Z per verificare se s è un quadrato perfetto è semplice, perchè basta considerare la somma modulo 2 degli esponenti, cioè essa deve essere pari. In Z[T] è un pò più complicato e dobbiamo svolgere diversi passi.

In pratica dobbiamo esaminare qualche altro teorema, che ci faciliti nel compito.

Teorema 7
Sia l in Z[T] tale da avere una fattorizzazione
\[l = (a_{1} + b_{1}T)^{e_{1}} (a_{2} + b_{2}T)^{e_{2}} ...\]
dove per ogni j, aj + bjT soddisfi il criterio di essere un fattore base algebrico. Se l è un perfetto quadrato in Z[T] allora \[\forall i,e_{i} \equiv 0 (mod 2).\]

In realtà l soddisferà anche qualche altra condizione utile.

Definizione simbolo di Legendre
Sia a in Z e p un numero primo. Il simbolo di Legendre sappiamo che è definito nel seguente modo:
\[\big(\frac{a}{b} \big) =\begin{cases} 1 \;  x^2 \equiv a mod p \\-1 \;  x^2 \not\equiv a mod p \\ 0 \;  p|a  \end{cases}\]

Teorema 8
Sia U un insieme finito di  coppie (a,b) tali che \[\prod_{{(a,b)} \in U} (a + bT)\] sia un quadrato perfetto in Z[T]. Allora per qualche coppia (s,q), con  q numero primo e s come nel Teorema 3 con \[(s,q)  \not\mid  a + bT\] per qualche (a,b) in U,
 \[\prod_{{(a,b)} \in U} \big(\frac{a+bs}{q} \big) = 1\]

 Da notare che dire:
 \[(s,q)  \not\mid  a + bT\]
 implica che:
 \[a \not\equiv -bs(mod q),\]
 così
 \[q \not\mid  a + bs\] e \[\big(\frac{a+bs}{q} \big) \neq 0\]

L'osservazione di sopra sarà utile successivamente.

I teoremi 7,8 danno una condizione necessaria ma non sufficiente per un elemento di Z[T] essere un quadrato perfetto.

Quindi per determinare se un l in Z[T] è un quadrato perfetto si lavora in questo modo:

a) per una fattorizzazione
 \[l = (a_{1} + b_{1}T)^{e_{1}} (a_{2} + b_{2}T)^{e_{2}} ... \]

 \[e_{j} \equiv 0 mod 2 \forall j.\]

b)  Sia Q un insieme di coppie di numeri (s,q) con q primo e s come da Teorema 3.

Scegliamo Q tale che \[(s,q) \not\mid a + bT\] per ogni a + bT che si verifica nella fattorizzazione di l.

 Verifichiamo che per ogni (s,q) in Q,
 \[\prod_{{(a,b)} \in U} \big(\frac{a+bs}{q} \big)= 1\] per U definito nel Teorema 8.

 L'insieme Q è chiamato carattere base quadratico e ogni (s,q) in Q è detto carattere quadratico.

 c) Se le precedenti due condizioni si verificano, allora l è probabilmente un quadrato perfetto in Z[T].

Nota che per incrementare questa probabilità, occorre incrementare il numero di elementi in Q. Nella seconda parte metteremo tutto assieme e poi faremo un esempio. Tosta questa spiegazione, e sicuramente se non arriviamo agli esempi vi sembrerà ancora il tutto poco chiaro.

Alla prox 

Riferimenti
[1]https://it.wikipedia.org/wiki/Polinomio_irriducibile
[2] Viet Pham Hoang - INTEGER FACTORIZATION WITH THE GENERAL NUMBER FIELD SIEVE
[3] Integer Factorization - Per Leslie Jensen - Master Thesis





giovedì 29 dicembre 2016

Seconda parte CFRAC in PARI GP

aggiornato al 27 gennaio 2017

E' il 29 dicembre, giornata da periodo natalizio. Da me fiocca la neve! Momento ottimale per rimanere un attimo in casa e fare un altro pezzetto di CFRAC.

Ho ritoccato un attimo il sorgente mostrato nel blog precedente e quindi riguardatelo un attimo. 

La seconda parte del software, cerca tutte le combinazioni (coppie, triple, etc) di righe della matrice E la cui somma modulo 2 per colonna dia un vettore con tutti zero (combinazioni lineari dipendenti).

L'algoritmo presentato non è ottimizzato e funziona bene per piccoli numeri n.

Ipotizzate che la parte che vi mostro è sempre contenuta e scritta in CFRAC.gp


\\ Trova le combinazioni semplici
\\ a gruppi di k
\\
\\ P è il vettore dei numeri pivot
\\ e ne contiene n_ele-1
\\ X è il vettore di nele elementi
\\ Y e Z sono di appoggio prima
\\ di inserire in Clist
\\ R. Turco


CombList(n_ele)=local();{
global(Clist);
Clist=listcreate(n_ele);
for(k=1,n_ele, \\ a k gruppi
   CombGroup(n_ele,k);
);
return(Clist);
}

CombGroup(n_ele,k)=local();{
global(Clist,Y,X,Z,P);
SetPivot(n_ele);

if(k==1,
   for(m=1,length(X),
       Y=listcreate(n_ele);
       listput(Y,X[m]); \\ ogni singolo elemento
       listput(Clist,Y);
   );
);
if(k==n_ele,
  listput(Clist,X); \\ tutti gli elementi in gruppo
);
if(k>1 && k!=n_ele,
    for(j=1,length(P),\\ per ogni elemento pivot
      Y=listcreate(n_ele);
      q = j;
      for(i=1,k-1,\\ predisponiamo i pivot
        if(q<=length(P),
           listput(Y,P[q]);
           q++;
        );
  );
       index=j+k-1;
       while(index<=length(X),\\ predisponiamo gli elementi non pivot
    Z=listcreate(n_ele);  
         for(f=1,length(Y),
           listput(Z,Y[f]);
         );
listput(Z,X[index]);
         listput(Clist,Z);
         index++;
      );

 \\ Aggiungo i pivot mancanti
 Y=listcreate(n_ele);
      q = j;
 v=k-1;
      for(i=1,k-1,\\ predisponiamo i pivot
        if(q<=length(P),
           listput(Y,P[q]);
  q=q+v;
  v=v+k+9;
        );
  );
  index=j+k;
  while(index<=length(X),\\ predisponiamo gli elementi non pivot
    Z=listcreate(n_ele);  
         for(f=1,length(Y),
           listput(Z,Y[f]);
         );
listput(Z,X[index]);
         listput(Clist,Z);
         index++;
      );
   );    
); \\ fine if

\\ pulizia della memoria
clean_memory();

\\ ordinamento di Clist
listsort(Clist,1);

}

SetPivot(n_ele)=local();{
  \\ P vettore dei pivot
   global(P,X);
   P=listcreate(n_ele-1);
   for(w=1,n_ele-1,listput(P,w));
   listsort(P,1);
   \\ X vettore di nele elementi
   X=listcreate(n_ele);
   for(w=1,n_ele,listput(X,w));
   listsort(X,1);
}

clean_memory()=local();{
  global(X,Y,P,Z);
  if(length(X)>0,listkill(X););
  if(length(Y)>0,listkill(Y););
  if(length(P)>0,listkill(P););
  \\if(length(Z)>0,listkill(Z););
}


SumMod2(E,n)=local(x=1,y=1);{

 CombList(matsize(E)[1]);
 v=vector(matsize(E)[2]); \\ vettore riferimento con zeri
 a=vector(matsize(E)[2]); \\ vettore delle somme modulo 2

 for(k=1,length(Clist), \\ Per tutte le combinazioni
   for(i=1,length(Clist[k]), \\ Prendo una combinazione
      a=(a+E[Clist[k][i],])%2;
      x = x*ListGcd[Clist[k][i]][2];
      y = y*ListGcd[Clist[k][i]][3];
 
   );
   if(a==v,
   \\ qua faccio il gcd
      r1= gcd(x - y,n);

     if(r1!=1 && r1!=n,
        print("factor = ", r1);
        print("factor = ", n/r1,"\n");
        return;
     );
     x=1;y=1;
   );
  );
  print("Sorry ...\n");
}

\\ CFRAC
\\ Output:
\\ Determina i fattori primi di un composto tabella
\\ Input:
\\ n è il numero da fattorizzare
\\ np è il numero di primi desiderati
\\ k è il k-esimo convergente della frazione continua

CFRAC(n,np,k)=local();{
if(isprime(n)==1, print("It's a prime number!\n"); return;);
 CreateMat(n,np,k);
 SumMod2(E,n);
}

Esempi di esecuzione

(12:33) gp > \r c:/CFRAC.gp

Esempio n=4141

(10:51) gp > CFRAC(4141,4,20)

k       pk%n    Bk+1    factor(B)
2       129     77      [7, 1; 11, 1]
3       193     -20     [-1, 1; 2, 2; 5, 1]
6       814     36      [2, 2; 3, 2]
8       3719    21      [3, 1; 7, 1]
11      2266    -84     [-1, 1; 2, 2; 3, 1; 7, 1]
12      3463    33      [3, 1; 11, 1]
13      232     -9      [-1, 1; 3, 2]
14      2570    5       Mat([5, 1])
15      2367    -84     [-1, 1; 2, 2; 3, 1; 7, 1]
17      3959    -4      [-1, 1; 2, 2]
18      3436    105     [3, 1; 5, 1; 7, 1]
19      3254    -21     [-1, 1; 3, 1; 7, 1]
20      3142    20      [2, 2; 5, 1]
factor base :
List([-1, 2, 3, 5, 7, 11])

List[k - pk%n - Bk+1]:
List([[2, 129, 77], [3, 193, 20], [6, 814, 36], [8, 3719, 21], [11, 2266, 84], [12, 3463, 33], [13, 232, 9], [14, 2570, 5], [15, 2367, 84], [17, 3959, 4], [18, 3436, 105], [19, 3254, 21], [20, 3142, 20]])
E.rows = 13

factor = 41
factor = 101

Esempio n=91
(10:51) gp > CFRAC(91,4,20)

k       pk%n    Bk+1    factor(B)
2       10      9       Mat([3, 2])
3       19      -3      [-1, 1; 3, 1]
5       33      -3      [-1, 1; 3, 1]
6       88      9       Mat([3, 2])
7       30      -10     [-1, 1; 2, 1; 5, 1]
8       27      1       matrix(0,2)
9       61      -10     [-1, 1; 2, 1; 5, 1]
10      88      9       Mat([3, 2])
11      58      -3      [-1, 1; 3, 1]
13      72      -3      [-1, 1; 3, 1]
14      10      9       Mat([3, 2])
15      82      -10     [-1, 1; 2, 1; 5, 1]
16      1       1       matrix(0,2)
17      9       -10     [-1, 1; 2, 1; 5, 1]
18      10      9       Mat([3, 2])
19      19      -3      [-1, 1; 3, 1]
factor base :
List([-1, 2, 3, 5, 11, 29])

List[k - pk%n - Bk+1]:
List([[2, 10, 9], [3, 19, 3], [5, 33, 3], [6, 88, 9], [7, 30, 10], [8, 27, 1], [9, 61, 10], [10, 88, 9], [11, 58, 3], [13, 72, 3], [14, 10, 9], [15, 82, 10], [16, 1, 1], [17, 9, 10], [18, 10, 9], [19, 19, 3]])
E.rows = 16

factor = 7
factor = 13

Esempio n=22703
(10:54) gp > CFRAC(22703,4,20)

k       pk%n    Bk+1    factor(B)
2       151     98      [2, 1; 7, 2]
7       8886    -38     [-1, 1; 2, 1; 19, 1]
9       3932    -119    [-1, 1; 7, 1; 17, 1]
11      2910    -119    [-1, 1; 7, 1; 17, 1]
13      4798    -38     [-1, 1; 2, 1; 19, 1]
18      3572    98      [2, 1; 7, 2]
20      14307   1       matrix(0,2)
factor base :
List([-1, 2, 7, 17, 19])

List[k - pk%n - Bk+1]:
List([[2, 151, 98], [7, 8886, 38], [9, 3932, 119], [11, 2910, 119], [13, 4798, 38], [18, 3572, 98], [20, 14307, 1]])
E.rows = 7

factor = 73
factor = 311

Esempio n=58807

Questo n non è un semiprimo, quindi dobbiamo fattorizzare in più fasi.
n = 7 * 8401 = 7 * 31 * 271

(13:34) gp > CFRAC(58807,4,20)

k       pk%n    Bk+1    factor(B)
2       243     242     [2, 1; 11, 2]
3       485     -3      [-1, 1; 3, 1]
4       19521   81      Mat([3, 4])
6       58804   9       Mat([3, 2])
7       39124   -27     [-1, 1; 3, 3]
9       57352   -27     [-1, 1; 3, 3]
10      52300   9       Mat([3, 2])
12      9       81      Mat([3, 4])
13      6561    -3      [-1, 1; 3, 1]
14      56611   242     [2, 1; 11, 2]
15      4365    -243    [-1, 1; 3, 5]
16      2169    1       matrix(0,2)
17      54442   -243    [-1, 1; 3, 5]
18      56611   242     [2, 1; 11, 2]
19      52246   -3      [-1, 1; 3, 1]
20      9       81      Mat([3, 4])
factor base :
List([-1, 2, 3, 11, 17])

List[k - pk%n - Bk+1]:
List([[2, 243, 242], [3, 485, 3], [4, 19521, 81], [6, 58804, 9], [7, 39124, 27], [9, 57352, 27], [10, 52300, 9], [12, 9, 81], [13, 6561, 3], [14, 56611, 242], [15, 4365, 243], [16, 2169, 1], [17, 54442, 243], [18, 56611, 242], [19, 52246, 3], [20, 9, 81]])
E.rows = 16

factor = 7
factor = 8401

(13:34) gp > CFRAC(8401,4,20)

k       pk%n    Bk+1    factor(B)
2       92      63      [3, 2; 7, 1]
4       275     16      Mat([2, 4])
5       2933    -135    [-1, 1; 3, 3; 5, 1]
7       4156    -120    [-1, 1; 2, 3; 3, 1; 5, 1]
9       1045    -105    [-1, 1; 3, 1; 5, 1; 7, 1]
10      8       64      Mat([2, 6])
13      6461    -48     [-1, 1; 2, 4; 3, 1]
14      1520    125     Mat([5, 3])
15      7981    -21     [-1, 1; 3, 1; 7, 1]
16      6981    160     [2, 5; 5, 1]
17      6561    -3      [-1, 1; 3, 1]
18      5794    40      [2, 3; 5, 1]
factor base :
List([-1, 2, 3, 5, 7])

List[k - pk%n - Bk+1]:
List([[2, 92, 63], [4, 275, 16], [5, 2933, 135], [7, 4156, 120], [9, 1045, 105], [10, 8, 64], [13, 6461, 48], [14, 1520, 125], [15, 7981, 21], [16, 6981, 160], [17, 6561, 3], [18, 5794, 40]])
E.rows = 12

factor = 31
factor = 271

Esempio n = 4780783 = 7^2 * 43 * 2269

Anche qui servono varie fasi e trovare k giusto:

(13:48) gp > CFRAC(4780783,12,10)

k       pk%n    Bk+1    factor(B)
3       4373    -3      [-1, 1; 3, 1]
4       1592865 729     Mat([3, 6])
5       3187915 -3638   [-1, 1; 2, 1; 17, 1; 107, 1]
6       4780780 9       Mat([3, 2])
7       3186460 -243    [-1, 1; 3, 5]
9       4767664 -27     [-1, 1; 3, 3]
10      4249828 81      Mat([3, 4])
factor base :
List([-1, 2, 3, 17, 23, 31, 47, 53, 67, 79, 101, 103, 107])

List[k - pk%n - Bk+1]:
List([[3, 4373, 3], [4, 1592865, 729], [5, 3187915, 3638], [6, 4780780, 9], [7, 3186460, 243], [9, 4767664, 27], [10, 4249828, 81]])
E.rows = 7

factor = 7
factor = 682969

(13:49) gp > CFRAC(682969,12,20)

k       pk%n    Bk+1    factor(B)
2       1653    533     [13, 1; 41, 1]
4       5785    744     [2, 3; 3, 1; 31, 1]
7       25619   -48     [-1, 1; 2, 4; 3, 1]
8       178160  1325    [5, 2; 53, 1]
9       203779  -297    [-1, 1; 3, 3; 11, 1]
11      514086  -720    [-1, 1; 2, 4; 3, 2; 5, 1]
12      141424  611     [13, 1; 47, 1]
13      113965  -248    [-1, 1; 2, 3; 31, 1]
15      540700  -923    [-1, 1; 13, 1; 71, 1]
16      682945  576     [2, 6; 3, 2]
18      85061   135     [3, 3; 5, 1]
20      87245   520     [2, 3; 5, 1; 13, 1]
factor base :
List([-1, 2, 3, 5, 11, 13, 31, 41, 47, 53, 61, 71, 73])

List[k - pk%n - Bk+1]:
List([[2, 1653, 533], [4, 5785, 744], [7, 25619, 48], [8, 178160, 1325], [9, 203779, 297], [11, 514086, 720], [12, 141424, 611], [13, 113965, 248], [15, 540700, 923], [16, 682945, 576], [18, 85061, 135], [20, 87245, 520]])
E.rows = 12

factor = 7
factor = 97567

(13:50) gp > CFRAC(97567,12,20)

k       pk%n    Bk+1    factor(B)
2       625     357     [3, 1; 7, 1; 17, 1]
3       937     -134    [-1, 1; 2, 1; 67, 1]
5       4373    -3      [-1, 1; 3, 1]
6       30544   282     [2, 1; 3, 1; 47, 1]
7       65461   -119    [-1, 1; 7, 1; 17, 1]
8       97254   402     [2, 1; 3, 1; 67, 1]
10      97564   9       Mat([3, 2])
11      64941   -94     [-1, 1; 2, 1; 47, 1]
12      96942   357     [3, 1; 7, 1; 17, 1]
13      64316   -243    [-1, 1; 3, 5]
14      63691   322     [2, 1; 7, 1; 23, 1]
17      84448   -27     [-1, 1; 3, 3]
18      31087   434     [2, 1; 7, 1; 31, 1]
20      67023   282     [2, 1; 3, 1; 47, 1]
factor base :
List([-1, 2, 3, 7, 17, 23, 31, 47, 53, 67, 79, 101, 103])

List[k - pk%n - Bk+1]:
List([[2, 625, 357], [3, 937, 134], [5, 4373, 3], [6, 30544, 282], [7, 65461, 119], [8, 97254, 402], [10, 97564, 9], [11, 64941, 94], [12, 96942, 357], [13, 64316, 243], [14, 63691, 322], [17, 84448, 27], [18, 31087, 434], [20, 67023, 282]])
E.rows = 14

factor = 43
factor = 2269

Avvertenza
L'algoritmo qui presentato è didattico e andrebbe ottimizzato ancora. Non sempre funziona per ogni n (probabilità 2/3 ricordate?) ed è molto lento per numeri molto grandi. Quando dà il messaggio "Sorry..." a volte è sufficiente aumentare leggermente, più volte, il valore np (per fB) o aumentare il valore k cioè il numero di convergenti ... Consiglio sempre di partire con np=4 e k=20, poi di provare ad alzare, più volte, il valore di np ed eventualmente k. Se non si fattorizza, allora il metodo non è del tutto adatto a quel numero. Nel frattempo però avrete saputo almeno se è un numero primo o meno.

Ad esempio con n=1162201201 occorre arrivare a np=15 e k=40 per iniziare ad avere un fattore 7 e proseguire in più fasi.

Più grande è n ed è meglio avere il factor base fB con numeri piccoli. Inoltre maggiore è n maggiore sarà la matrice E da ottenere e il metodo di riduzione di Gauss potrebbe essere pesante, mentre il metodo a blocchi di Lanczos, salvando su file parti di matrice, potrebbe essere una soluzione migliore da farsi in un linguaggio come il C o il Fortran.

In pratica si parte da "relazioni parziali" che magari non portano a combinazioni linearmente dipendenti, mentre allargando il fattore base fB si riesce a ottenere una "relazione totale" che ci porta rapidamente al quadrato e quindi alla fattorizzazione. Non sempre si riesce facilmente a fattorizzare. Provate ad esempio 387400807. Se ci riuscite avvisatemi. Si sa che è 387400807 = 19441· 19927 ma provateci con CFRAC.

Brrr .... fa freddo!
Musica classica di Mozart, un bel caminetto, un whisky ed un toscanello sono degli ottimi antitodi!

Nei prossimo post puntiamo oltre, cioè al General Number Field Sieve. E qui bisognerà lavorare un poco per comprendere bene il metodo.

p.s. Nel frattempo Buon 2017!!!



Alla prox



martedì 27 dicembre 2016

Prima parte di CFRAC in PARI GP

aggiornato al 5 gennaio 2017

Mi avete scritto di essere interessati a vedere l'algoritmo CFRAC automatizzato. Qui vi propongo la parte  fino alla soluzione rappresentata dalla matrice E degli esponenti dopo l'eliminazione di Gauss modulo 2.


Listato algoritmo CFRAC prima parte
\\ Output:
\\ Determina un factor base fB per CFRAC
\\ Input:
\\ n è il numero da fattorizzare
\\ np è il numero di primi desiderati
\\

fb(n,np)=local(j=1);{

  global(fB);

  fB=listcreate(np);
  \\ aggiungo due valori di default
  listput(fB,-1);
  listput(fB,2);

  \\ cerco gli altri primi che soddisfano il simbolo di Legendre (residui quadratici)
  forprime(i=2,n,
    if(kronecker(n,i)==1 && j<=np,
 listput(fB,i);
 j++;
);
  );
  \\ ordino la lista
  listsort(fB,1);

  return(fB);
}

\\ListSearch Cerca un elemento nella lista
\\ ritorna la posizione se la trova o 0
ListSearch(ele,mylist)=local(i=1);{
   while(i<=length(mylist),
     if(ele==mylist[i],
   return(i);
);
i++;
   );
   return(0);
}

\\ Output:
\\ Determina la tabella di elementi da considerare rispetto a fB CFRAC
\\ Input:
\\ n è il numero da fattorizzare
\\ np è il numero di primi desiderati
\\ k è il k-esimo convergente della frazione continua
\\ In output si ottiene stampa della tabella e restituzione lista dei fattori della tabella

GMod2(n,np,k)=local(no=1);{

  \\ global comodo per fare test
  global(ListFactorTab,ListGcd);

  mylistfb=fb(n,np);
  a=contfrac(sqrt(n));
  ck=contfracpnqn(a,k);
 
  ListFactorTab = listcreate(k);
  ListGcd = listcreate(k);
  print("\nk","\t", "pk%n","\t", "Bk+1","\t", "factor(B)");
  for(i=1,length(ck),
      B=(ck[1,i]^2-n*ck[2,i]^2);
      m=factor(B);
 for(j=1,matsize(m)[1],
    if(ListSearch(m[j,1],mylistfb)==0,
   no=0;
);
 );
 if(no!=0,
   print(i,"\t", ck[1,i]%n,"\t", B,"\t", factor(B));
listput(ListFactorTab,factor(B));
listput(ListGcd,[i,ck[1,i]%n,abs(B)]);
 );
 no=1;
  );
  print("factor base :");
  print(fB);
  print("\nList[k - pk%n - Bk+1]:\n",ListGcd);
  return(ListFactorTab);
}

\\ Output:
\\ Crea la matrice E soluzione da AE=b
\\ con metodo di eliminazione Gauss modulo 2
\\ sugli esponenti
\\ a partire dalla tabella ListfactorTab
\\ Input:
\\ n è il numero da fattorizzare
\\ np è il numero di primi desiderati
\\ k è il k-esimo convergente della frazione continua

CreateMat(n,np,k)=local(size=0,j=1, i=1,pos=0);{

  global(E);

  GMod2(n,np,k);
  size=length(ListFactorTab);
  E=matrix(size,np+2);
  print("E.rows = ",size,"\n");
  for(k=1,size,
    for(j=1,matsize(ListFactorTab[k])[1],
     pos=ListSearch(ListFactorTab[k][j,1],fB);
if(pos>0,
   E[i,pos]=ListFactorTab[k][j,2]%2;
pos=0;
);
  );
  i++;
  );
  return(E);
}

Esempio di esecuzione
Supponiamo che il sorgente si chiami CFRAC.gp. Useremo il numero composto n=4141, vogliamo altri 4=np primi (per un totale di 6) nel fattore base (fB) e 20 convergenti e riuscire ad osservare la tabella risultante.

Ok via!

gp> \r c:/CFRAC.gp

gp > CreateMat(4141,4,20)
2       129     77      [7, 1; 11, 1]
3       193     -20     [-1, 1; 2, 2; 5, 1]
6       814     36      [2, 2; 3, 2]
8       3719    21      [3, 1; 7, 1]
11      2266    -84     [-1, 1; 2, 2; 3, 1; 7, 1]
12      3463    33      [3, 1; 11, 1]
13      232     -9      [-1, 1; 3, 2]
14      2570    5       Mat([5, 1])
15      2367    -84     [-1, 1; 2, 2; 3, 1; 7, 1]
17      3959    -4      [-1, 1; 2, 2]
18      3436    105     [3, 1; 5, 1; 7, 1]
19      3254    -21     [-1, 1; 3, 1; 7, 1]
20      3142    20      [2, 2; 5, 1]
E.rows = 13

%32 =
[0 0 0 0 1 1]

[1 0 0 1 0 0]

[0 0 0 0 0 0]

[0 0 1 0 1 0]

[1 0 1 0 1 0]

[0 0 1 0 0 1]

[1 0 0 0 0 0]

[0 0 0 1 0 0]

[1 0 1 0 1 0]

[1 0 0 0 0 0]

[0 0 1 1 1 0]

[1 0 1 0 1 0]
[+++]
(16:22) gp >

Se avrò tempo vi mostrerò anche la seconda parte ...

ps. il simbolo [+++] indicato da PARI GP vi avverte che ce ne sono altre di righe e che per default fa vedere solo 12 righe.

Se ad esempio fate E[13,1] oppure ListFactorTab[13][1,1] osserverete che c'è un'altra riga,
mentre con ListFactorTab vedrete la lista di tutti i fattori invece.

Infatti Listsize(ListFactorTab)
dà 13

Buon Natale e alla prox!







sabato 17 dicembre 2016

CFRAC - Quadrati congruenti e frazioni continue con PARI GP



I numeri RSA per la crittografia , quelli proposti su Wikipedia, sono stati scelti accuratamente e non possono essere fattorizzati con gli algoritmi classici di Trial Division, Fermat, Pollard, ECM.

Questi numeri RSA possono essere fattorizzati più rapidamente con degli algoritmi della famiglia "Number Field Sieve". Il termine "più rapidamente" in realtà è relativo e dipende sempre sia dal numero, sia da situazioni di parallelismo/potenza computers e comunque si tratta solitamente di mesi.

Per arrivare al GNFS (General Number Field Sieve) inizieremo a spiegare, qui e in altri post, diverse cose, semplificando e dando per scontato qualche Teorema, a vantaggio del filo conduttore. Tuttavia quando avrete compreso il quadro generale, dovrete consolidare le parti date per scontate.

Non è semplice semplificarvi tutto. Ci provo.

Nel seguito indicherò con a|b per dire che a è divisore di b. Mentre con il simbolo != intenderemo diverso.

Congruenze Quadrate
Il tutto ha avuto inizio da un'idea di Fermat, rivisitata da Legendre.

Se è soddisfatta la congruenza x^2 = y^2 mod n con 0<=x<=y<=n e x!=y e x+y!=n allora vuol dire che:

n | (x-y)(x+y)  e se n=pq allora:

pq | (x-y)(x+y) => p|(x-y) o p|(x+y)  => q|(x-y) o q|(x+y) l'or qui è esclusivo

Facendo gcd(n,x-y) o gcd(n,x+y) si può ottenere un fattore non banale e, quindi fattorizzare n.

Si possono avere varie combinazioni possibili che riassumeremo nella tabella successiva.

p|(x-y)     p|(x+y)    q|(x-y)     q|(x+y)    gcd(n,x-y)       gcd(n,x+y)      fattore non banale
si no si no n 0 no
si no            no            si p q si
si no si si n q si
no si si no q p si
no            si            no si 0 n no
no si si si            q                    n                    si
si si            si no n p                    si
si si no si p n                    si
si si            si            si            n                    n                    no


Abbiamo in tabella 9 "casi totali" con 6 "si (successo)", quindi abbiamo 2/3 di probabilità di trovare un fattore primo non banale con questo metodo.

Tutti gli algoritmi che vedremo hanno i seguenti step:
- trovare un insieme di relazioni che sono B-smooth sopra un fattore base
- risolvere il sistema di equazioni lineari e cercare le relazioni con i quadrati trovati
- calcolare il gcd tra il composto n e i quadrati trovati

Numeri lisci o B-smooth
B-Smooth (numero B-liscio) che vuol dire? Se n è un intero positivo, allora n è un numero B-smooth se tutti i suoi fattori primi in cui esso si scompone, presi con esponente unitario,  {p1,p2,p3,...,pmax} sono tali che ogni pi<=B, dove B=pmax.

Frazioni continue
Legendre iniziò a valutare l'idea delle farzioni continue per la fattorizzazione ma solo con Morrison e Brillhart si ebbe la definizione definitiva del metodo CFRAC.

CFRAC era in voga negli anni '70-'85 e usato per i numeri di Fermat del tipo Fn. Con esso fu fattorizzato nel '70 il numero \[ F_{7}=2^{2^{7}} + 1 = 340282366920938463463374607431768211457 =\]
\[59649589127497217  *   5704689200685129054721 \]

CFRAC è di interesse ancora nel progetto ProthSearch.

Una frazione continua è del tipo:
\[x_{k}=[a_{0},a_{1},a_{2},a_{3},...]\]
\[x_{k} = a_{0} + \frac{b_{1}}{a_{1}+\frac{b_{2}}{a_{3}+..}}\]

Se i vari bi=1 allora si tratta di frazioni continue semplici.

Ora xk è detto "convergente k-esimo" di x se ci siamo fermati al k-esimo termine.

Ad esempio se x=3.51 e usate contfrac di PARI GP esce x=[3,1,1,24,2]. Qui i termini sono solo k=5.

Una frazione continua infinita
\[[a_{1},...,a_{k},,...]\]

non è altro che il limite del convergente \[x_{k}=[a_{1},...,a_{k}]\]

Teorema
Un numero reale x si può rappresentare sempre come frazione continua.

Algoritmo di Euclide
Per ottenere una frazione continua ci si basa sull'algoritmo di Euclide e la frazione continua si ottiene come nel seguito:

\[x_{0} = x\]
\[q_{0} = floor(x_{0})\]
\[x_{1} = \frac{1}{x_{0}-q_{0}}\]
\[q_{1} = floor(x_{1})\]
\[x_{2} = \frac{1}{x_{1}-q_{1}}\]
\[. . . . . . ........\]
\[q_{i} = floor(x_{i})\]
\[x_{i+1} = \frac{1}{x_{i}-q_{i}}\]

Teorema 
Se consideriamo x=sqrt(n), con n free square cioè n non quadrato, allora x è espandibile in una frazione continua periodica:
 \[x=\sqrt{n} = [a0,\overline{a1,...,an,2a0}]\]

Definizione
I convergenti Pn/Qn di una frazione continua semplice [q0,q1,q2,q3,...] sono dati da:

\[\frac{P_{0}}{Q_{0}} = \frac{q_{0}}{1}
\frac{P_{1}}{Q_{1}} = \frac{q_{0}q_{1} +1}{q_{1}}
. . .
\frac{P_{i}}{Q_{i}} = \frac{q_{i}P_{i-1} + P_{i-2}}{q_{i}Q_{i-1} + Q_{i}-2} ,i>=2\]

CFRAC usa proprio i convergenti della radice di n.

I convergenti che sono B-smooth su un fattore base fB vengono presi come relazioni.

Le relazioni sono messe in uno step di algebra lineare e quelle relazioni che producono un quadrato perfetto, i fattori del convergente sono posti nel gcd (Massimo comun divisore). Se il gcd con n è diverso da 1 e n si ottiene un fattore primo non banale e, quindi la fattorizzazione.

Riassumiamo l'algoritmo e proviamo a fare un esempio.

Algoritmo CFRAC
input : numero composto n.
output : un fattore p non banale di n.

Step 1: (Cercare i convergenti)
Calcolare l'espansione di frazione continua della radice di n e ottenere le relazioni (Pk, Bk+1) e il fattore base fB.

Step 2: (Cercare i quadrati delle relazioni con i fattori dei convergenti con metodi di algebra lineare)
Sulle relazioni usiamo il metodo di eliminazione di Gauss (o quello di Strassen più efficiente), in modo da risolvere il sistema di equazioni lineari.

Step 3: (gcd)
p può essere cercato con p = gcd(n,x±y). Abbiamo 2/3 di probabilità di riuscirci.

Se non troviamo nulla allora si ricomincia dallo Step 1, aumentando il numero di elementi di fB e/o variando il numero di convergenti.

Il Metodo
Per lo Step 1, i convergenti li possiamo ottenete con contfracpnqn(a,k) dove a=contfrac(n) e k è un valore da noi scelto, per avere il k-esimo convergente.

Usate help con ?? contfracpnqn e vedete gli esempi.

I convergenti Pk/Qk sono forniti da contfracpnqn come due vettori: quello di sopra è Pk, quello di sotto è Qk.

Per lo Step 2 occorre aver chiaro alcune cose.

Se si prende Pk/Qk di radice di n allora è:

Pk^2  - n*Qk^2 = (-1)^(k+1) * Bk+1     (a)

che implica che:
Pk^2 = (-1)^(k+1) * Bk+1 mod n (b) 

Per ottenere le relazioni, l'idea è prendere coppie (Pk, Bk+1) per produrre quadrati.

Ricordiamo che si ottiene un quadrato perfetto se le potenze dei fattori sono tutte pari.

Per cercare i prodotti di Bk che danno quadrati perfetti facciamo prima la loro fattorizzazione e poi li combiniamo così che gli esponenti diventano pari.

La fattorizzazione, poichè sono piccoli numeri, la possiamo fare col Trial Division.

Selezioniamo un insieme di primi dai fattori di ogni Bk.

Se p|Bk allora Pk = n*Qk^2 mod p

In tal caso n è residuo quadratico modulo p.

Ricordiamo che il simbolo di Legendre (n/p) vale 1 se n è un residuo quadratico modulo p, vale -1 se non lo è, vale 0 se p divide n.

Per cui selezioneremo un insieme di primi q tale che il simbolo di legendre (n/q) = 1. Questo insieme di primi è il fattore base fB.  Ci dovremmo sicuramente mettere il -1, per rendere positivi i risultati, e il 2 per iniziare a considerare l'elenco dei primi consecutivi.

Esempio 

n=4141;
a=contfrac(sqrt(n))
[64, 2, 1, 5, 1, 3, 3, 3, 5, 1, 4, 1, 3, 14, 25, 1, 2, 31, 1, 5, 6, 3, 1, 2, 1, 4, 2, 2, 2, 2, 4, 1, 2, 1, 3, 6, 6]
ck=contfracpnqn(a,20)
[64 129 193 1094 1287 4955 16152 53411 283207 336618 1629679 1966297 7528570 107366277 2691685495 2799051772 8289789039 259782511981 268072301020 1600144017081 9868936403506]

[1 2 3 17 20 77 251 830 4401 5231 25325 30556 116993 1668458 41828443 43496901 128822245 4036986496 4165808741 24866030201 153361989947]

La prima riga sono i Pk. Dobbiamo trovarci i Bk. Usiamo la (a) con PARI GP. I fattori dati da PARI GP leggiamoli come una lista separata da punto e virgola e dove il primo termine è un numero primo e il secondo l'esponente.

Usiamo il comando PARI GP:
for(k=1,20,print(k,"\t", ck[1,k]%n,"\t", B=(ck[1,k]^2-n*ck[2,k]^2),"\t", factor(B)));

Ecco la tabella:
k Pk      Bk+1 Factors
1       64      -45      [-1, 1; 3, 2; 5, 1]
2       129     77      [7, 1; 11, 1]
3       193    -20      [-1, 1; 2, 2; 5, 1]
4       1094   87      [3, 1; 29, 1]
5       1287  -31      [-1, 1; 31, 1]
6       814     36       [2, 2; 3, 2]
7       3729  -37      [-1, 1; 37, 1]
8       3719   21      [3, 1; 7, 1]
9       1619  -92      [-1, 1; 2, 2; 23, 1]
10      1197   23     Mat([23, 1])
11      2266  -84     [-1, 1; 2, 2; 3, 1; 7, 1]
12      3463   33     [3, 1; 11, 1]
13      232     -9      [-1, 1; 3, 2]
14      2570    5      Mat([5, 1])
15      2367  -84     [-1, 1; 2, 2; 3, 1; 7, 1]
16      796     43     Mat([43, 1])
17      3959   -4      [-1, 1; 2, 2]
18      3436   105   [3, 1; 5, 1; 7, 1]
19      3254   -21    [-1, 1; 3, 1; 7, 1]
20      3142    20     [2, 2; 5, 1]

Dobbiamo scegliere ora un fattore base fB. Quale scegliamo?

Il metodo è:

a) provare i primi q tali che il simbolo di Legendre (n/q)=1, per trovare i residui quadratici,ì. Useremo la funzione PARI GP Kronecker(n,q). Il simbolo di Legendre e quello di Jacobi sono una generalizzazione di quello di Kronecker.

b) individuare solo i primi q desiderati. Troppi primi rallentano il metodo.

Per non fare i calcoli a mano (di solito i numeri n sono enormi), per trovare fB candidato possiamo usare un programmino scritto ad hoc del tipo:

\\ Determina un factor base fB (output)
\\ input:
\\ n è il numero da fattorizzare
\\ k è il piccolo numero di primi desiderato


fb(n,k)=local(j=1);{
  global(fB);
  listkill(fB);
  fB=listcreate(k);
  listput(fB,-1);
  listput(fB,2);
  forprime(i=2,n,
    if(kronecker(n,i)==1 && j<=k,
 listput(fB,i);
 j++;
);
  );
  listsort(fB,1);
  print(fB);
}

eseguiamo:
\r fb(4141,4)
risultato del factor base (sono 4 perchè abbiamo per default aggiunto -1 e 2):
fB=[-1, 2, 3, 5, 7, 11]

Abbiamo scelto k=4 per avere al massimo 6 termini. Maggiore è il numero di termini scelto, più termini ci sono nel sistema delle equazioni da risolvere e questo si traduce in una maggiore lentezza nella ricerca dei fattori non banali.

Il -1 ci servirà per compensare il segno dei Bk negativi. Il 2 perchè ci serve per avere i pari.

A questo punto la tabella di prima va rivista, per escludere i Bk+1 che non contengono elementi di fB=[-1, 2, 3, 5, 7, 11], ad esempio quelli per k=4,5,7,9,10,16.

La tabella rimanente è:
k Pk     Bk+1 Factors
2       129     77      [7, 1; 11, 1]
3       193     -20     [-1, 1; 2, 2; 5, 1]
6       814     36      [2, 2; 3, 2]
8       3719    21     [3, 1; 7, 1]
11      2266    -84   [-1, 1; 2, 2; 3, 1; 7, 1]
12      3463    33    [3, 1; 11, 1]
13      232     -9      [-1, 1; 3, 2]
14      2570    5      Mat([5, 1])
15      2367    -84   [-1, 1; 2, 2; 3, 1; 7, 1]
17      3959    -4     [-1, 1; 2, 2]
18      3436    105  [3, 1; 5, 1; 7, 1]
19      3254    -21   [-1, 1; 3, 1; 7, 1]
20      3142    20    [2, 2; 5, 1]


Ora siamo al passo di algebra lineare. Vediamo cosa fare.

Supponiamo che la fattorizzazione dei Bk è del tipo:
\[B_{1} = p_{1}^{a_{11}}*p_{2}^{a_{12}}*...p_{r}^{a_{1r}}
...
B_{s} = p_{1}^{a_{s1}}*p{2}^{a_{s2}}*...p_{r}^{a_{sr}}\]

dove p1=-1 e p1,...,ps sono i primi del fattore base.

Noi abbiamo bisogno di trovare degli esponenti e1,...,es che valgano 0 o 1 tali che:
B1^e1*B2^e2...*Bs^es   sia un quadrato perfetto.

Da qui deve essere:

B1^e1*B2^e2...*Bs^es  = (p1^a11 * p2^a12 * ... * pr^a1r)^e1 ...\(p1^as1 * p2^as2 *... * pr^asr)^er

Quindi ci serve trovare degli e1,...,es che per ogni indice i sia pari l'espressione

e1*a1i+e2*a2i+ ... + es*asi

Cioè

a11 e1} + a12 e2 +...+a1s es = 0 mod 2

.....

ak1 e_{1} + ak2 e2 +...+aks es = 0 mod 2

In algebra lineare corrisponderebbe a:

A * e = 0 mod 2  (3)

dove A è la matrice degli aij mentre e=(e1,e2,...,es).

In PARI GP la soluzione e dalla (3) si può ottenre con matsolvemod

Sulla matrice A di dimensione nxm (può essere n=m oppure no) applichiamo il metodo di eliminazione di Gauss modulo 2.

Creiamoci adesso la matrice E col metodo di eliminazione di Gauss modulo 2, guardando se gli esponenti di A sono pari. I vettori riga vi della matrice E li segniamo rispetto ai primi di fB:

Matrice E delle soluzioni
-1 2 3 5 7 11
 0 0 0 0 1 1             v1
 1 0 0 1 0 0 v2
 0 0 0 0 0 0             v3
 0 0 1 0 1 0             v4
 1 0 1 0 1 0             v5
 0 0 1 0 0 1            v6
 1 0 0 0 0 0             v7
 0 0 0 1 0 0 v8
 1 0 1 0 1 0             v9
 1 0 0 0 0 0             v10
 0 0 1 1 1 0             v11
 1 0 1 0 1 0             v12
 0 0 0 1 0 0             v13

Le soluzioni del sistema di equazioni produrranno delle combinazioni

pi1^2 pi2^2 ... pik^2 = Bi1+1  Bi2+1  ... Bik+1 mod n

dove Bi1+1  Bi2+1  ... Bik+1 è un quadrato.

Dalla matrice E ottenuta in sostanza ci interessano le combinazioni linearmente dipendenti. Ne abbiamo diverse da provare e cioè quelle la cui somma degli 1 in colonna è pari:

v2+v4+v5+v13

.....
v5 + v12

Ad esempio
 v2+v4+v5+v13
 1 0 0 1 0 0
 0 0 1 0 1 0
 1 0 1 0 1 0
 0 0 0 1 0 0

Non è detto che immediatamente la soluzione dia l'esito sperato ma con vari tentativi si riesce a ottenere (ad esempio non con la prima combinazione lineare).

Ora dobbiamo combinare i Pk trovati nella tabella sottraendo dei fattori tra quelli di fB

Con la prima combinazione lineare è:
v2+v4+v5+v13
gcd(193·3719·2266·3142-20*3*2*7,n)=4141    KO

Ad esempio con vari tentativi si ottiene con
v5+v12
1 0 1 0 1 0
1  0 1 0 1 0

gcd(2266*3254-2*3*7,n)=41
gcd(2266*3254+2*3*7,n)=101

Infatti 4141=41*101

Ovviamente l'esempio è stato fatto  per un n piccolo.

Ora potreste anche automatizzare il metodo con un programmino CFRAC :-)
Alla prox

Riferimenti
http://www.ams.org/journals/mcom/1975-29-129/S0025-5718-1975-0371800-5/S0025-5718-1975-0371800-5.pdf