/*********************************************************
Languasco & Zaccagnini, Introduzione alla Crittografia, Hoepli.
Capitolo 3.8.4 pag. 72 e
Capitolo 9.6.3 pag. 278
*********************************************************/

{ MillerRabin(n,c=1) = local(m,t,s,a,b,e,c1);
if((n <= 3)&(n>1), print(n, " e` primo"); return);
if( n<=1, error("Input non valido"));
if(!bittest(n,0), print(n, " e` pari"); return(0));
/*********************************************************
 INIZIALIZZAZIONE: Posto t = n-1, s = 0, e fino a 
               che t e` pari si pone t = t/2 e s = s+1 
*********************************************************/
m = n-1;
t = m;
s = 0;  
while(!bittest(t,0), t=t >>1; s++);         
/*********************************************************
     alla fine del while si ha  n-1=2^s * t, t dispari 
*********************************************************/ 
c1=c;
while(c > 0,
/*********************************************************
TEST: usando un generatore di numeri pseudocasuali,
          si sceglie un a tale che 1 < a < n. Si 
          pone e = 0, b = a^t (mod n). Se b == 1 o 
          b == n-1 si pone c = c-1 e  si ripete il test. 
*********************************************************/
     until(lift(a) > 1,  a = Mod(random(n),n));
     e = 0; 
     b = Mod(lift(a),n)^t;    
     if(b == 1 | b==m, c--; next);        
/*********************************************************
   n e` pseudoprimo forte in base a;  ripartenza del ciclo  
*********************************************************/
      while( b != 1 && b != m && e <= s-2,  
                          b = Mod(lift(b),n)^2; e++);
      if(b != m, print(n, " e` composto");return(0)); 
/*********************************************************
QUADRATI :  Se b != (+-)1 (mod n) e e <= s-2, 
          calcoliamo b = b^2 (mod n) e e = e + 1. 
          Se b != n - 1, n e` composto, ritorno 0 e 
          l'algoritmo termina, altrimenti n e` 
          pseudoprimo forte in base a, decremento 
          di c e ripartenza del ciclo   
*********************************************************/ 
c--);    
print(n," e` pseudoprimo forte in ", c1, " basi a (random)");
return(1);
}