Salq(minimo=1,maximo=10) = local();{
 
 for(i=minimo,maximo,
     Alq(i);
 );

}


Alq(n,op="+",m=0,k=1) = local(ret=1, str="", str2="");{

 str = concat(str,n);
 str = concat(str,",");

 i = 0;
 retp = n;

 while(ret!= 0,  

       i++;
       vec = divisors(n);
       l = length(vec);
       retp2 = ret;
       ret = sum(i=1, l-1, vec[i] );
       

       if( op = "+", ret = ret*k + m; );
       if( op = "-", ret = ret*k - m; );
       if( op = "*", ret = ret*k + m; );
       if( op = "^", ret = ret^k + m; );

       if( retp == ret,     
           ret = 0;
       );
       if( retp2 == ret,
           retp2 = -1;  
           if( retp2==-1, str=concat(str,ret); str=concat(str," cycle"););  
           ret = 0;
           
       );

       if( retp != ret & retp2 != -1,    
           if(ret==0,   str= concat(str,ret);); 
           if(ret!=0,   str= concat(str,ret); str=concat(str,",");); 
           n = ret;
       );
       
 
       if( ret==0 & retp2 != -1, str = concat(str," - len = "); str=concat(str,i));
       if( ret==0 & i==1 & retp2 != -1, str = concat(str," - perfect number"););
       if( ret==0 & i==2 & l > 2 & retp2 != -1, str = concat(str," - amicable number");); /* Not prime number */
       if( ret==0 & i>2 & retp2 != -1, str = concat(str," - sociable number"););
       
       
 );
 print(str);

}

AlqSum(n)=local(ret=1);{

 vec = divisors(n);
 ret = sum(i=1, length(vec)-1, vec[i] );
 return(ret);

}

/*
*  Verifica se un numero pari è perfetto
*/
IsPrf(a) = local(ret=0);{
 
 if( a%2 != 0, error("Perfect numbers are even!"); return(0););
 if( n == AlqSum(n), 
     print(a, " is a perfect number");
     ret=1;
 );
 if( n != AlqSum(n), 
     print(a, " isn't a perfect number");
     ret=0;
 );

 return(ret);
}

/*
*  Cerca num perfetti dispari
*/
SPrf(n1=1, n2=1000)=local();{

 for(i=n1,n2,

     if( i%2 == 0, i++ );
     s=AlqSum(i);
     if( s == i & i%2 != 0, print(s, " Num. perf. dispari!!!"););
     num=i;
 );
 print("Numeri esaminati ", num-1 );
 return;
}

/*
*
*  Cerca coppie amicabili con tecnica Al-Faris
*/

Samic(n1=2, n2=1000)=local();{

 if( n1<2, error("deve essere n1>=2"););
 for(k=n1,n2,

     p=(3*(2^(k-1)))-1;
     q=(3*(2^(k)))-1;
     r=(9*(2^(2*k-1)))-1;
     if(isprime(p) & isprime(q) & isprime(r) == 1,
         N=(2^k)*(p*q);
         M=(2^k)*r;
         print("(",N,",",M,")");
     );
 );   
 print("Numeri esaminati ", n2-n1+1 );
 return;
}

/*
*
*  Cerca tutte le coppie amicabili 
*/

Allamic(n1=2, n2=1000)=local(s=0, prec1=0);{

 if( n1<2, error("deve essere n1>=2"););

 for(k=n1,n2,
     s=AlqSum(k);
     if( AlqSum(s) == k && k! = s & k != prec1,
         print("(",k,",",s,")");
         prec1=s; 
    );
     
 );   
 print("Numeri esaminati ", n2-n1+1 );
 return;
}