/*
**  R. Turco 2010
**  versione 3.1
**
*/


/*

** Simple Method
*/


pMod(base,exponent,modulus)= local(ret=1); { 
    
ret = lift(Mod(base, modulus)^exponent);
    
return(ret);

}



/*

**  Fast Mod Exp: exponentiation by squaring 

**
**  It uses "Bruce Schneier's Algorithm" 

**  (Applied Cryptography, 2e, ISBN 0471117099)

**
**  Concrete implementation by R. Turco

*/


ModExp(base,exponent, modulus) = local(ret=1,i=1); {
   
 b = binary(exponent);
   
 i= matsize(b)[2];
 while(exponent>0 & i>0,
         
        if(b[i]==1,
             
           ret = lift(Mod((ret * base), modulus));
     	 
        );	 
		 
        exponent = exponent >> 1;
	     
        base = lift(Mod((base * base), modulus));
		 
        i--;
    
 );
    
 return(ret);

}





/*
** Square & Multiply 
**  
*/
SqMod(base,exponent,modulus)= local(residue=1); { 
    
 b = binary(exponent);
  
 len= matsize(b)[2];
 for(i=1,len,
      residue=(lift(Mod(residue,modulus)^2)*(base^b[i])) % modulus;
 );
 return(residue); 
}





/*
*   Primality Fermat-Carmichael
*   
*
*/
FC(number, tec=1)= local(ret=0); { 
    
 ret = witness(number, 40);
 if( ret == 1,
     if( tec ==1, ret= pMod(2,number-1, number););
     if( tec ==2, ret= SqMod(2,number-1, number););
     if( tec ==3, ret= ModExp(2,number-1, number););
     if( ret != 1, ret=0;);
 );
 return(ret); 

}

/*
*  Witness test
*
*/

witness(n,nbasi)=local(ret=0,exponent=0,t=0, u=0, status=2); { 
 
 
 exponent = n-1;
 while(exponent%2==0,
       exponent = exponent/2;
       t = t+1;
 );
 u = (n-1)/(2^t);
 q=2;
 while( q < nbasi & q < n & status!=0, 
        val = q;
        q++;
        exponent = u;
        i = t;
        while(i>-1,
          ret=ModExp(val,exponent,n);
          if( ret==1, status=1; break;);
          i--; 
          val=ret;
          exponent=2;
   
        );
        if( ret!=1, status=0; break;);
 );
 if( status == 2, status = 0; );
 return(status);
}


/*
*  Repunit Primality
*
*  2, 19, 23, 317, 1031, 49081, 86453, 109297
*/

Repunit(p) = local(ru=0); { 
   ru=(10^p-1)/9;
   return(ru);
}

/* 
   it returns: 
   0 if x is composite
   1 if x is prime
 
   Test Di Maria Giovanni  
        
*/

LRFT(x) = local(ret=0, rep=0, res=0); { 
   trap( , return(0),
     res=lift(Mod(2,x)^(eulerphi(x)/(x*2)))%1000;   
   );
   if( res==0, ret=1);
   return(ret);
}

/* 
   it returns: 
   0 if x is composite
   1 if x is prime
 
   Test Rosario Turco   
        
*/
LRFT2(x) = local(ret=0, r=0, e=0); {
   r = Repunit(x);
   e = (r-1)/x;
   trap(,return(0),
         ret=ModExp(3,e,r)%1000;
         ret2 = ret%1000;
   );
   if(ret2 == 0, 
      ret=1;
   );
   if(ret2!=0, 
      ret=0;
   );
   return(ret);
}



lambda(n)= local(i=0,j=0); { 
   len=matsize(factor(n));
   f=matrix(len[1],len[2]);
   f=factor(n);
   v=vector(len[1]);
   if( len[1] == 1 & f[1,1]==2,
     v[1] = (f[1,1])^(f[1,2]-2);
   );
   if(f[1,1]!=2,
    for(j=1,len[1],
          v[j] = (f[j,1]-1)*(f[j,1]^(f[j,2]-1));
       );
   );
   print("factors : ",f);
   print("lambda  : ",lcm(v));

   return(lcm(v));
 }


lambda1(n)= local(i=0,j=0); { 
   len=matsize(factor(n));
   f=matrix(len[1],len[2]);
   f=factor(n);
   v=vector(len[1]);
   for(j=1,len[1],
     v[j] = f[j,1]-1;
   );
   return(lcm(v));
 }

ex(inf,sup)=local(i=0); {
   for(i=inf,sup, 
       a=Repunit(i);
       b=lambda1(a);
       c=factor(a);
       write1("k:\\work\\repunit2.txt"," digit: ", i, " Repunit: ", a, " lambda1: ", b, " resto: ", a%b, " factor: ", c, "\n");
   );
    
 }