/*
*    LibThN
*
*    Library "Number Theory"
*
*    Author: 	Rosario Turco 
*    Copyright:	(C) 2007-2020
*    Version: 	4.4
*    Date:	02/09/09
*
*    Product: 	PARI/GP  version 2.4.2
*
*    Blog:	http://MATHBuildingBlock.blogspot.com
*
*    Site:	www.geocities.com/SiliconValley/Port/3264
*    Papers: 	Database CNR Solar
*
*
*   
*    SOFTWARE WITHOUT ANY WARRANTY WHATSOEVER.
*
*
*/

/********************************************************
* GetHelp
* It lists the functions
* 
* 
********************************************************/

{GetHelp() = local(); 
  print1("   List of helps\n");
  print1(" # HelpThN()		- help on the Theory of Numbers\n");
  print1(" # HelpFactOrCrypt()	- help on the Factorization or Crypting\n");
  print1(" # HelpCollatz()	- help on the Collatz's Problem\n");
  print1(" # HelpMath()		- help on the Math\n");


}

{HelpThN() = local();
  print1("   HELP on Theory of Numbers\n");
  print1(" # listLandauPrimes(n) [It writes on file in c: the Landau’s primes up to n]\n");
  print1(" # getLandauPrimes(n) [It counts the Landau’s primes up to n]\n");
  print1(" # getLandauPrimesFrom(n,np,L) [It counts from np to n, with L=0 or previous calculation L] \n"); 
  print1(" # listTwinPrimes(n) [It writes on file in c: the twin primes up to n]\n");    
  print1(" # getTwinPrimes(n) [It counts the twin primes up to n]\n");    
  print1(" # listPolignacPrimes(n,d) [It writes on file in c: the Polignac's primes up to n with d distance]\n"); 
  print1(" # getPolignacPrimes(n,d) [It counts the Polignac's primes up to n with d distance]\n"); 
  print1(" # listSophieGermainPrimes(n) [It writes on file in c: the Sophie Germain's primes on file in c:]\n");
  print1(" # getSophieGermainPrimes(n) [It counts the Sophie Germain's primes up to n ]\n");  
  print1(" # listGold(n) [It writes on the screen the pair of solution of Goldbach] \n");
  print1(" # getGold(n) [It counts the number of pair of solutions of Goldbach] \n");
  print1(" # getLevy(n) [It counts the Levy’s pair primes such that for N>5 is N = p +2q where p and q are prime numbers] \n");
  print1(" # listLevy(n) [It printss the Levy’s pair primes such that for N>5 is N = p +2q where p and q are prime numbers] \n");
  print1(" # getMatLevy(n) [It returns a Levy's Matrix of numeric progressions p+2q, where p and q are prime numbers and n is the nth prime number] \n");
  print1(" # getMsngLevy(n) [It returns missing odd numbers in Levy's Matrix - n is the nth prime number] \n");
  print1(" # getDetLevy(n) [It returns the Determinant in Levy's Matrix] \n");
  print1(" # listSigma(n1,n2) [It returns sigma(n) (sum of divisors) from n1 to n2] \n");
  print1(" # listIntRobin(n1,n2) [It returns the integers which sigma(n) < exp(gamma)*n*log(log(n)) from n1 to n2, n1>1 (Robin's criterion)] \n");
  print1(" # listRatioRobin(n1,n2) [It returns r=(exp(gamma) * n * ln(ln(n)))/sigma(n) from n1 to n2]\n");
}

{HelpFactOrCrypt() = local();
  print1("   HELP on Factorization\n");
  print1(" # getRSA(1,n)  [It creates a composite that is the product of two primes, n is the numebr of digits of one prime]\n");
  print1(" # crackRSA(n) [n factors or 2 factors: with n factors you must repeat on the single factor]\n");
  print1(" # scompFactor(n) [n factors]\n");
  print1(" # cryptWordRSA(Myword): In RSA MyWord must be a string in quotes\n");
  print1(" # decryptWordRSA(N, puK, secret): N, puK and secret are output of cryptWordRSA\n");
  print1(" # cryptWordElG((mc, p, a, gb): EL GAMAL-> View the source\n");
  print1(" # decryptWordElG(secretElG, p, b, ga): EL GAMAL-> View the source\n");
   
 print1(" # qsfRSA(n) [2 factors: it must be a semiprime RSA\n");
}

{HelpCollatz() = local();
  print1("   HELP on 3x+1\n");
  print1(" # G(n) [It writes on the screen the Collatz's sequence] \n");
  print1(" # G1(n) [It returns the numbers of step in the Collatz's sequence]\n");
  print1(" # findMr(n1=0,n2=0, r=0.0, W=0, Y=0) [It finds Max r=Pos/G(N) in (n1,n2)]\n");
  print1(" # Gr(n) [It writes on the screen the Collatz's (3x+1)/2] \n");
  print1(" # G1r(n)  [It returns the numbers of step in the Collatz's (3x+1)/2]\n");
  print1(" # isopm(n) [It returns the "isopath max"]\n");
  print1(" # isog(n,n1,n2) [It searches G(n)=G(x) for x from n1 to n2]\n");
  print1(" # gminus(n,n1,n2) [It searches G(n1)=n1-n from n1 to n2]\n");
  print1(" # gplus(n,n1,n2) [It searches G(n1)=n1+n from n1 to n2]\n");
  print1(" # gbiuniv(n1,n2) [It searches G(a)=b e G(b)=a]\n");
  print1(" # detectCnum(n) [It returns the type of the number]\n");
  print1(" # f2k(n) [It returns the number as sum of power of 2]\n");
  print1(" # maxk2k(n) [It returns 2^k less than n or equal to n]\n");
  print1(" # isnCollatz(n) [It returns 1 if n is a Collatz's number (2^k-1)/3]\n");
  print1(" # isnBizzarre(n) [It returns 1 if n is a Bizzarre's number (2^k-4)/4]\n");
  print1(" # Vp(n) [It returns the parity vector of n]\n");
  print1(" # P(n,k) [It returns the number of Parity of order k]\n");
  print1(" # D(n,k) [It returns the number of 1 in the parity vector]\n");
  print1(" # T(n,k) [It returns the approssimation of Sk value of order k]\n"); 
  print1(" # maxOrd(n) [It returns the max order of the parity vector]\n");

}

{HelpMath() = local();
  print1("   HELP on Math\n");
  print1(" # log2(n) [n: it must be a integer\n");
  print1(" # logb(n,b) [n: it must be a integer, b: base (integer)\n");
  print1(" # loge(n) [natural log - n: it must be a integer\n");

 }

/*********************************************************
* Classic Collatz's Problem (3x+1)
*
* Rosario Turco
*
*********************************************************/
{listCollatz(n) = local (j, k, f, s="");

 if( n<2 , error("listCollatz(n): You must insert an integer n>1");); 

 print1("\n", n, " ");

 j=0;
 k=0;
 f=0;
 M=0;
 A=0;

 while( n>1,
        f=0;
        if( n%2 != 0 & f==0, n=3*n+1; f=1;);
        if( n%2 == 0 & f==0, n=n/2; f=1;);
        if( n > M, M=n; A=j+2;);
        j++; 
        print1(n, " ");
      );

 s=concat(s,"\n\nG(N): ");
 s=concat(s,j);
 s=concat(s, " Max: ");
 s=concat(s,M);
 s=concat(s, " Pos: ");
 s=concat(s,A);
 s=concat(s," r=Pos/G(N): ");
 a= (1. * A)/j;
 s=concat(s, a);
 print(s); 
}

{getCollatz(n) = local (j);

 if( n<2 , error("getCollatz(n): You must insert an integer n>1");); 
 j=0;
 while( n>1,
        f=0;
        if( n%2 != 0 & f==0, n=3*n+1; f=1;);
        if( n%2 == 0 & f==0, n=n/2; f=1;);
        j++; 
      );
 return(j);
}

{G(n)  = listCollatz(n);}
{G1(n)  = getCollatz(n);}

/*
** findMr(n1,n2)
**
** It finds Max r=Pos/G(N) in (n1,n2)
** in input r, W e Y can be useful to remember previous values for a new interval
**
**
*/

{findMr(n1=0,n2=0, r=0.0, W=0, Y=0)= local (j, k, f, M, A, s="");

 if( n1<5 | n2 < 5 | n2 <=n1, error("findG(n1,n2): You must insert n1>4, n2>4 and n2>n1");); 

 for(q=n1,n2, 
 j=0;
 k=0;
 f=0;
 M=0;
 A=0;
 n = q;

 while( n>1,
        f=0;
        if( n%2 != 0 & f==0, n=3*n+1; f=1;);
        if( n%2 == 0 & f==0, n=n/2; f=1;);
        if( n > M, M=n; A=j+2;);
        j++; 
      );
 
 
 print1("\nN=", q);
 s=concat(s,"\nG(N): ");
 s=concat(s,j);
 s=concat(s, "\t Max: ");
 s=concat(s,M);
 s=concat(s, "\t Pos: ");
 s=concat(s,A);
 s=concat(s,"\t r=P/G(N): ");
 a= (1. * A)/j;
 if( r < a, r=a; W=q; Y=M;);
 s=concat(s, a);
 print(s);
 s=""; 
 );
 print("\n rmax = ", r, " N: ", W, " Max: ", Y);
 print("\n");
}


/*
** List of numbers of the Glide 
** with 3x+1 algorithm restated (by Terras)
**
** Gr(N)
** G1r(N)
**
*/
{Gr(n) = local (j, k, f, s="", p=0, i=0);

 if( n<2 , error("Gr(n): You must insert an integer n>1");); 

 s=concat(s,"S0=");
 s=concat(s,n);
 s=concat(s, " ");

 j=0;
 k=0;
 f=0;
 p=n; /* start value */

 while( n>1,
        f=0;
        if( n%2 != 0 & f==0, n=(3*n+1)/2; f=1;);
        if( n%2 == 0 & f==0, n=n/2; f=1;);
        if( n < p, i=j;);
        s=concat(s,"S");
        s=concat(s,j+1);
        s=concat(s,"=");
        s=concat(s,n);
        s=concat(s, " ");
        j++; 
      );
 print(s);
 print("glide: ", j+1);
 print("stopping time: ", i, "\n");
 
}

{G1r(n) = local (j, k, f);

 if( n<2 , error("G1r(n): You must insert an integer n>1");); 

 j=0;
 k=0;
 f=0;

 while( n>1,
        f=0;
        if( n%2 != 0 & f==0, n=(3*n+1)/2; f=1;);
        if( n%2 == 0 & f==0, n=n/2; f=1;);
        j++; 
      );
 return(j+1);
}


/*
**  parity Vector
**
*/
{Vp(n,t=0) = local (i=0, j=1, m=0, f, v);

 if( n<2 , error("Vp(n): You must insert an integer n>1");); 

 i=G1r(n);
 v=vector(i);
 v[j]=n%2;
 if( t==1, s=concat("V", j-1); s=concat(s,"="); s=concat(s,v[j]);print(s););
 j++;

 while( i>j-1,
        f=0;
        if( n%2 != 0 & f==0, n=(3*n+1)/2; f=1;);
        if( n%2 == 0 & f==0, n=n/2; f=1;);
        v[j]=n%2;
        if( t==1, s=concat("V", j-1); s=concat(s,"="); s=concat(s,v[j]);print(s););
        j++; 
      );

 return(v); 
}

{maxOrd(n) = local();
  return(maxk2k(n)+2);
}

/*  P(n, k) 
**  Parity Number of n of order k (as a binary number of parity vector)
**  
**
*/

{P(n, k=-1, t=1) = local (i=0, j=1, p=0, v);

 if( n<2 , error("P(n,k): You must insert an integer n>1");); 

 i=G1r(n);

 if( k<0 , k=i-1;);  /* K is order from 0 to i-1 */

 if( k>i-1 , error("P(n,k): You must insert an integer k valid"));

 
 v=vector(i);

 v=Vp(n);

 while( j<k+2,
        p = p + v[j]*(2^(j-1));
        if( t==1, s=concat("P", j-1); s=concat(s,"="); s=concat(s,p);print(s););
        j++; 
      );

 return(p); 
}

/*
**  D(n,k) , k default value is the max order
**
*/
{D(n, k=-1) = local (i=0, j=1, p=0, v);

 if( n<2 , error("D(n,k): You must insert an integer n>1");); 

 i=maxOrd(n);

 if( k<0 , k=i-1;);  /* K is order from 0 to i-1 */

 if( k>i-1 , error("D(n,k): You must insert an integer k valid"));

 v=vector(i);
 v=Vp(n);
 
 while( j<k,
        p = p + v[j];
        j++; 
      );

 return(p); 
}


/*
**  T(n,k), k default value is the max order.
**
*/
{T(n,k=-1) = local (p=0.0, i=0);

 if( n<2 , error("T(n,k): You must insert an integer n>1");); 
 i=maxOrd(n);
 if( k<0 , k=i-1;);  /* K is order from 0 to i-1 */

 if( k>i-1 , error("T(n,k): You must insert an integer k valid"));

 p = 1. * n * 3^D(n,k) * 2^(-k);
 return(p); 
}


/********************************************************
*
* getLevy(N)  
* It counts the Levy’s pair primes such that 
* for N>5 is N= p +2q where p and q are prime numbers
*
* listLevy(n)
* It prints the Levy’s pair primes such that 
* for N>5 is N= p +2q where p and q are prime numbers
*
********************************************************/

{getLevy(n) = local(p,j);
  
  if( Mod(n,2) == 0, error("You must insert odd numbers")); 
  if( n<6, error("You must insert n>5")); 

  j=0;
  for(p=2,n,   
     if( isprime(p), 
         if( n - p > 2,   
             if( Mod(n-p,2) == 0, 
                 if(isprime((n-p)/2),j=j+1);
             );
         );
     ); 
  );
  return(j);
}
{listLevy(n) = local(p);
  
  if( Mod(n,2) == 0, error("You must insert odd numbers")); 
  if( n<6, error("You must insert n>5")); 


  for(p=2,n,   
     if( isprime(p), 
         if( n-p > 2 , 
             if( Mod(n-p,2) == 0, 
                if(isprime((n-p)/2), print1(p); print1("+"); print1(n-p); print1("  "););
             );
         );
     ); 
  );

}


/********************************************************
* listLandauPrimes(n)    
* It returns the Landau’s primes up to n on a file
* in directory c:
*
* getLandauPrimes(n)  
* It counts the Landau’s primes up to n
*
********************************************************/
{listLandauPrimes(n) = local(i, j);

  j=0;
  for(i=1,n,   
     if( isprime(i*i + 1), write1("c:\\landau.txt"," ", i*i + 1));
     if( isprime(i*i + 1), j=j+1); 
     if( i*i+1 > n, break;); 
  );
  print1(" #Landau's prime up to ", n, " : ", j); 

}
{getLandauPrimes(n) = local(i,j);

  j=0;
  for(i=1,n,   
     if( isprime(i*i + 1), j=j+1); 
     if( i*i+1 > n, break;); 
  );
  return(j);
}
{getLandauPrimesFrom(n,np,L) = local(i,j, a);

  j=L;
  a=floor(np);
  for(i=a, n,   
     if( isprime(i*i + 1), j=j+1); 
     if( i*i+1 > n, break;); 
  );
  return(j);
}


/********************************************************
* listTwinPrimes(n)    
*
* It returns the Twin Primes up to n on a file
* in directory c:
*
* getTwinPrimes(n)  
* It counts the Twin Primes up to n
*
* 
********************************************************/
{listTwinPrimes(n) = local(i, j, k);

  j=0;
  for(i=3,n-2,
      
     k = i + 2;

     if( isprime(i), 
         if( isprime(k), write1("c:\\twinprimes.txt"," ", i, "-", k); j=j+1);
             i = i + 2;
     );  
     if( k > n, break;); 

   );
  print1(" #Twin Primes up to ", n, " : ", j); 

}

{getTwinPrimes(n) = local(i,j,k);
  j=0;
  for(i=3,n-2,
      
     k = i + 2;

     if( isprime(i), 
         if( isprime(k),j=j+1);
             i = i + 2;
     );  
     if( k > n, break;); 
   );
  return(j);
}

/********************************************************
* listPolignacPrimes(n,d)    
* It returns the Polignac's Primes up to n with distance d
*
* getPolignacPrimes(n, d)  
* It counts the Polignac's Primes up to n with distance d
*
* 
********************************************************/
{listPolignacPrimes(n,d) = local(i, j, k);

  if( d==0, print("help: GetPolignacPrimes(n,d)"); return);
  if( Mod(d,2) != 0, error("You must insert even numbers")); 

  j=0;
  for(i=3,n-d,
      
     k = i + d;

     if( isprime(i), 
         if( isprime(k), write1("c:\\polignacprimes.txt"," ", i, "-", k); j=j+1);
             i = i + d;
     );  
     if( k > n, break;); 

   );
  print1(" #Polignac Primes  up to ", n, " : ", j); 

}

{getPolignacPrimes(n,d) = local(i,j,k);

  if( d==0, print("help: GetPolignacPrimes(n,d)"); return);
  if( Mod(d,2) != 0, error("You must insert even numbers")); 
  
  j=0;
  for(i=3,n-d,
     k = i + d;

     if( isprime(i), 
         if( isprime(k),j=j+1);
             i = i + d;
     );  
     if( k > n, break;); 
   );
  return(j);   
}

/********************************************************
* listSophieGermainPrimes(n,d)    
* It returns the Sophie Germain's primes up to n 
*
* getSophieGermainPrimes(n)  
* It counts the Sophie Germain's primes up to n
*
* 
********************************************************/
{listSophieGermainPrimes(n) = local(i, j, s);

  j=0;
  for(i=2,n,
     if( isprime(i), 
       s = 2*i + 1; 
       if( isprime(s), write1("c:\\sophiegermain.txt"," ", i,":" s);j=j+1;);
     );  
     if( i > n, break;); 
   );
  print1(" #Sophie Germain Primes  up to ", n, " : ", j); 
}
{getSophieGermainPrimes(n) = local(i, j, s);

  j=0;
  for(i=2,n,
     if( isprime(i), 
       s = 2*i + 1; 
       if( isprime(s), 
           j=j+1;
         );
     );  
     if( i > n, break;); 
   );

   return(j);


}

/*********************************************************
* Goldbach
*
° It returns the pairs (x,n-x) of prime numbers, 
° not repeated, equals to an even number n.
* 
*********************************************************/


{listGold(n) = local ();

  if( n<=4, 
      error("listGold(n): n>4 and even");
  );

  if( Mod(n,2) != 0, 
      error("listGold(n): You must insert an even n>4");
    ); 

  for(p=3,n, 
      if(isprime(p) & isprime(n-p) & p <= n-p, 
         print1("(", p, "," , n-p, ")", " ");
      );
   );
       
}


{getGold(n) = local(s=0,p=0); 

  if( n<=4, 
      error("getGold(n): n>4 and even");
  );  

  if( Mod(n,2) != 0, 
      error("getGold(n): You must insert an even n>4");
  ); 
  
  s=0;

  for(p=3, n, 
      if(isprime(p) & isprime(n-p) & p <= n-p, 
         s++;
      );
  );

  return(s);
}

/*
*
*   It returns all solutions of Goldbach in a matrix 2xn
*
*/

{memGold(n) = local (mat,col, j);

  if( n<=4, 
      error("memGold(n): n>4 and even");
  );

  if( Mod(n,2) != 0, 
      error("memGold(n): You must insert an even n>4");
  ); 

  col = getGold(n);
 
  mat = matrix(2, col);

  j=1; 
  for(p=3,n, 
      if(isprime(p) & isprime(n-p) & p <= n-p, 
         mat[1,j] = p; 
         mat[2,j] = n - p; 
         j++;
      );
  );
  return(mat);     
}

/*
*   It chooses, in way random, a solution between all 
*   solutions of Goldbach and inserts it in a vector
*
*/

{randGold(n) = local (mat,col, j, vec);

  if( n<=4, error("randGold(n) --> with n>4 and even"););

  if( Mod(n,2) != 0, error("randGold(n) --> You must insert an even n>4");); 

  col = getGold(n);

  mat = matrix(2, col);

  mat = memGold(n);

  j = random(col);

  if( j==0, j++);
  if( j>col, j--);
 
  vec = vector(2);

  vec[1]=mat[1,j];
  vec[2]=mat[2,j];

  return(vec);     
}


/*
*   It chooses the last solution of Goldbach and inserts it in a vector
*
*/

{lastGold(n) = local (mat,col, j, vec);

  if( n<=4, error("lastGold(n) --> with n>4 and even"););

  if( Mod(n,12) != 0, error("lastGold(n) --> You must insert an even number N=12n");); 

  col = getGold(n);

  mat = matrix(2, col);

  mat = memGold(n);

  j = col;
 
  vec = vector(2);

  vec[1]=mat[1,j];
  vec[2]=mat[2,j];

  return(vec);     
}
/*
*    It generates two prime numbers, with Goldbach, 
*    but it rejects N=12n so it doesn't use 
*    twin prime numbers 
*    because RSA could be easily cracked
*
*    The dimension of a random value determines 
*    the number of digits.
* 
*
*/
{getPrimes(n) = local (j, vec, test);

  vec=vector(2);
  test=1;

  j=random(10^n); /* numbers of digit that you want */

  while( test>0 | j<5 , j++; 
         if( !(j%2==0) | (j%12==0), j++;);
         if( j>4 & (j%2==0) & !(j%12==0), test=0; break; );
        );  

  vec=randGold(j);

  return(vec);     
}


/*
*    It generates two prime numbers, without Goldbach. 
*
*    The dimension of a random value determines 
*    the number of digits.
*    n is the number of digits 
*
*/
{getPrimesFast(n) = local (p, q, vec);

  vec=vector(2);
  p = nextprime(random(10^n)); 

  q = nextprime(random(10^(n+5))); 
  vec[1]=p;
  vec[2]=q;
  return(vec);     
}

/*
*    getRSA  
*    It generates a semiprime RSA with Goldbach 
*
*    If we use getRSA(1,n) 
*    with 1 we obtain a trace of two prime numbers in a semiprime
*    n is the number of digit that you want with one prime of the product 
*    the compost will have between 2*n-1  to 2*n digits
*
*/

{getRSA(f,n,al=0) = local (vec, valore);

  vec=vector(2);
  if( al == 0,   vec=getPrimesFast(n);); /* without Goldbach */   
  if( alg == 1, vec=getPrimes(n););      /* with Goldbach */
  
  valore=vec[1]*vec[2];
  if( f==1, print("\n getRSA: p = ", vec[1], " q = ", vec[2], " p*q = ", valore););
  return(valore);     
}

/*
* Factorization Area
*
*/

/*
*   crackRSA(n)
*
*   It prints two prime factors of a sempiprime RSA
*   such that n=p*q
*
*   It uses few lines of code 
*
*   Use:
*
*   crackRSA(getRSA(1))
*   or 
*   crackRSA(n)
*   
*   It uses a strategy MPQS
*
*
*   Note: If you use with getRSA(1,n) we raccomend 
*   n max 25 digit, otherwise you obtain the message
*   floor: precision too low in truncr (precision loss in truncation).
*   in crackRSA
*       
*   This is a didactic software.
*   Function factor in PARI is faster.
*
*/

{crackRSA(p) = local (s, vec, qp);

  vec=vector(2);

  qp = 4*p;

  s=floor(sqrt(qp)); 
  
  while( !(s^2>=qp) | (issquare(s^2 - qp) == 0) | (frac(s^2 - qp)!= 0.0), 
         s++; 
       ); 

  vec[1]=floor((s+sqrt(s^2-qp))/2);
  vec[2]=floor((s-sqrt(s^2-qp))/2);

  return(vec);     
}

/*
*  If it returns [0], is a prime number.
*  scompFactor uses bigomega 
*
*  This is a didactic software.
*  Function factor in PARI is faster.
*  
*/
{scompFactor(p) = local ( col, f1, f2, f3, s, i, jOut, vecIn, vecOut, vecStack);

  vecIn=vector(2);

  col = bigomega(p);

  vecOut=vector(col);
  vecStack=vector(2*col);

  print("\nscompFactor");
  print("\n#fattori previsti: ", col); 
 
  f1=0;
  f2=0;
  f3=0;
  
  s=1;
  jOut=1;
  iStack=1;

  while( s<col, 

         vecIn=crackRSA(p);

         if( isprime(vecIn[2]), vecOut[jOut]=vecIn[2]; jOut++; f1=1; print("->",vecIn[2]););  
         if( isprime(vecIn[1]), vecOut[jOut]=vecIn[1]; jOut++; f2=1; print("->",vecIn[1]););  

        
         if( s<col & f2==1 & f1==0, p=vecIn[2];);
         if( s<col & f1==1 & f2==0, p=vecIn[1];);     
          
         if( s<col & f1==0 & f2==0, vecStack[iStack]=vecIn[1]; 
            p=vecIn[2]; iStack++; vecIn[1]=0; vecIn[2]=0; f3++;); 

         if( iStack>1 & s>1 & f1==1 & f2==1 & (vecStack[iStack-1]!= 0), 
             p=vecStack[iStack-1]; f3--; vecStack[iStack-1]=0; iStack--; 
         );

         if( iStack==1 & s>1 & f1==1 & f2==1 & (vecStack[iStack]!= 0), p=vecStack[iStack]; f3--; 
             vecStack[iStack]=0;
         );

           
         if(p==1, s=col;);

         f1=0;
        f2=0;
         s++; 
         
       ); 

  return(vecOut);     
}

/*
*  Cryptography Area
*
*/

/************************************************************
* Am example of Crypt/Decrypt in RSA with public key e private key
*
* Out alphabet is of 26 char:
* A, B, C, D, E, F, G, H, I, J, K, L,M, N, O, P, Q, R, S, T, U, V,W, X, Y, Z
*
* Example of Word :
* HELLOWORLDBYROSARIOTURCO
*
* Attention: max digits of text you can crypt with private key must be: log n / log 26
* Text of 20 chars for time
*
* Example of use in sequence:

  createCryptRSA()                       [here -> KPr=(N, prk) KPu=(N,puk)] 
  get26Ascii("HELLOWORLDROSARIOTURCO")   [here -> m ]
  cryptWordRSA(m, N, prK)                [here -> secret] 
  decryptWordRSA(N, puK, secret)         

  or
  testRSA("HELLOWORLD")
*  
*
* 
* Rosario Turco
*************************************************************/


{createCryptRSA() = local (p,q, KP, phi);
 
 KP=vector(3);

 p = nextprime(random(10^40)); 

 q = nextprime(random(10^35)); 
 

 KP[1]=p*q;             /* N */
 
 phi = (p-1) * (q-1);  /* phi N */

 print("phi : ", phi);


 /* 
 *  Person A produce own private key and public key */  
 /* Public key = (N, puk) */
 /* Private key = (N, prk) */

 /* I put them in a vector KP for convenience in the test*/

 puK = random(q);

 print("\nN   : ", KP[1]);
 
 while(gcd(phi,puK)!=1,puK=puK+1);

 KP[2]=puK;  
 print("puK : ", puK);


 prK = lift(Mod(puK,phi)^(-1));


 KP[3]=prK;
 print("prK : ", prK);
 
 
 /* I verify that prk*puk is congrue con 1 mod phi */
 test = (prK * puK) % phi;

 if( test == 0, 
     error("Test NOK: It isn't congrue 1 mod phi"); );

 print(" ");

 return(KP);

}

/*
*  External Function
*/

{get26Ascii(str) = local (MyW, weight, i, a, m);

 
 weight = length(str);

 MyW=Vecsmall(str);  


 m=0;  

 /* 
    I use a power of 27 otherwise char Z will have 
    26*26^weight and this doesn't work well
    with getMsg 
 */

 for( i=1,weight,

      a = MyW[i]-65+1;
      m = a*(27^weight) + m;
      weight--; 
  );

 
 return(m);
 
}


/*
*   External function
*/
{getMsg(num) = local (a=num,d=27,r=0, msgM="", amax=0);

 print("num :", num);

 print("d :", d);
 while( d>26,

        amax = maxk26k(a);
        d = floor(a/(27^amax));
        r = a%(27^amax);
        msgM=concat(msgM,Strchr(d+64));

        if( r>=26, a = r; d=27;);
        if( r<26 & r!=0, msg=concat(msgM,Strchr(r+64)); d=1;);
        
      ); 

  return(msgM);
 
}

/* 
*  Internal function
*/

{maxk26k(n) = local (k);

 k=0;

 while( 27^k <= n, k++;);

 return(k-1);

}

/*
*   External function
*/

{cryptWordRSA(m, N, puK) = local (secret);
 
 
 secret = lift(Mod(m,N)^puK); 
 print("\nsecret value: ", secret);
 print(" ");
 
 
 return(secret);
 

}
/*
*   External function
*/

{decryptWordRSA(secret, N, prK) = local (desecret, msg);
 
 desecret =  lift(Mod(secret,N)^prK);

 print("decrypt value: ", desecret);
  
 print ("\nTest OK if decrypt value is equal to m value");
 
 return(desecret);


}

/*
   Test function 

   Person B transmit a message to A with public key di A 
   Person A decrypt with own private key

   If person A transmit a message to B must use the public key of B



*/

{testRSA(strMyWord) = local (vec,ms,secretRSA,desecretRSA,msg);
 
 /* B has got a KeyPub and transmit to A  a message */
 Key=vector(2);
 Key=createCryptRSA();

 ms = get26Ascii(strMyWord);
 print("m   : ", ms); 

 secretRSA=cryptWordRSA(ms, Key[1], Key[2]);


 /* Now A decrypt */
 
 desecretRSA=decryptWordRSA(secretRSA,Key[1],Key[3]);

 msg = getMsg(desecretRSA);
 print("\nmsg decriptato: ", msg);


 return;


}


/* 
*  EL GAMAL
*
*  External function
*
*/

{createCryptElG() = local (KPEG, p, g, a, b, ga, gb);
 
 KPEG=vector(6);  /* It's Public Key of El Gamal */

 p = nextprime(random(10^40)); 

 g = random(p);

 a = random(p);  /* It's Private Key of El Gamal user A */                

 b = random(p);  /* It's Private Key of El Gamal user B */                

 ga=lift(Mod(g,p)^a);
 gb=lift(Mod(g,p)^b);
 

 /* I put them in a vector KPEG for convenience in the test*/
 
 KPEG[1]=g;
 KPEG[2]=a;   /* a for user A and b for user B   */ 
 KPEG[3]=ga;  /* ga for user A and gb for user B */
 KPEG[4]=p;   /* p is shared */
 KPEG[5]=b;
 KPEG[6]=gb;


 return(KPEG);

}

/*
*   External function
*/

{cryptWordElG(mc, p, a, gb) = local (msecret);
 
 
 msecret = lift(mc*Mod(gb,p)^a); 
 print("\nsecret value: ", msecret);
 print(" ");
 
 
 return(msecret);
 

}

/*
*   External function
*/

{decryptWordElG(secretElG, p, b, ga) = local (desecretElG, msg);

 
 desecretElG =  lift(secretElG*(Mod(ga,p)^b)^-(1));

 print("decrypt value: ", desecretElG);
  
 return(desecretElG);


}

{testElG(strMyWordTest) = local (msc,Msecret,Mdesecret,msgs);
 
 /* B has got a KeyPub and transmit to A  a message */
 Key=vector(6);

 Key=createCryptElG();

 /* B tx to A */
 msc = get26Ascii(strMyWordTest);
 print("msc   : ", msc); 

 
 Msecret=cryptWordElG(msc, Key[4], Key[2],Key[6]);

 /* A decrypts */
 Mdesecret=decryptWordElG(Msecret,Key[4],Key[5],Key[3]);

 print("Decodifica\n");
 msgs = getMsg(Mdesecret);
 print("\nmsg decriptato: ", msgs);
 

 return;


}


/*
********************************************************
* Quadratic Sieve Factoring 
*	
*     x^2 = a^2 mod N
*     x^2 = a^2 + k*N (k=0, k=1)
* 
* Variant of Fermat's algorithm
* This method is good for odd numbers
*
* Rosario Turco
********************************************************
*/

{qsfRSA(n) = local (col, vecOut, xq, a, x1, x2);

 if( n<2 | frac(n) != 0.0 , error("qsfRSA(n): You must insert an integer n>1")); 

 col = bigomega(n);
 
 /* trivial solution: 1 and n it is a prime*/
 if( col == 1, print1("qsfRSA: It's a prime number: ", n); return(););

 vecOut=vector(col);

 if( col != 2, print1("qsfRSA: It isn't a RSA number. You must insert a RSA number"); return();); 

 /* if it is a even number ... */
 if( Mod(n,2)==0, vecOut[1]=2; vecOut[2]=floor(n/2); return(vecOut););

 xq=0;  /* xq = x^2 */
 a=0;

 while( (!ispower(xq,2) | xq == 0), a++; xq=a^2+n;
 );

 x1=a;
 x2=floor(sqrt(a^2+n));

 vecOut[1]=x2-x1; 
 vecOut[2]=x2+x1;

 /* if n is a perfect square */
 if( x2-x1==1 & x2+x1==n, vecOut[1]=floor(sqrt(n)); vecOut[2]=vecOut[1];);

 return(vecOut);

}


/*********************************************************
* Shor Factoring 
*	
*     n
*     f(x) = a^x mod n
* 
* Variant of "Shor's algorithm"
*
* Rosario Turco
*
* RSA
*          0         0     5   
* 39460021111111111107165109=2207*55234		r=
* 1178888888888888888771=1439*39095		r=
* 11788888811=		   			r=
* 2467757=2417*1021				r=
* 691217=1021*677				r=57460
* 99097=41*2417 				r=12080
* 4171=43*97    				r=336
*********************************************************/
{shorRSA(n, flag) = local (a, vecOut, i, r, j, initseq, val1, val2, repeat);

 if( n<2 | frac(n) != 0.0 , error("shorRSA(n): You must insert an integer n>1")); 

 /*
 if( bigomega(n) == 1, print1("shorRSA: It's a prime number: ", n); return(););
 if( bigomega(n) > 2, print1("shorRSA: You must insert a RSA number!"); return(););

 */
 vecOut=vector(2);

 /* trivial case */
 if( Mod(n,2)==0, vecOut[1]=2; vecOut[2]=floor(n/2); return(vecOut);); 
 if( Mod(n,3)==0, vecOut[1]=3; vecOut[2]=floor(n/3); return(vecOut);); 
 if( ispower(n,2)==1, vecOut[1]=floor(sqrt(n)); vecOut[2]=vecOut[1];return(vecOut); );
 
 /* Strategy to choose the value a */
 repeat = 1;
 a=2; 

 while( repeat == 1 & a<n,
 /* if r odd then repeat = 1 */
 while( repeat == 1,

 repeat = 0;
 r = 0;

 /* We find a<n and gcd(a,n)=1 */
 while((gcd(a,n)!=1), a++;);
  

  /* The first of the sequence */
 initseq = a%n;  /* a^1 mod n */
 
 i=2;


 /* We search the lenght of sequence jstep from the second element i=2 */
 while( initseq != (a^i)%n, r++; if(flag==1, print1("."); i++; ));
 
 r++;

 if(r%2 != 0, repeat=1; print1("d", r););


 ); 

 j=r/2;

 val1=(a^j)-1;
 val2=(a^j)+1;

 if(flag == 1, print1("\n\nr: ", r); print1("\nj: ", j););

 vecOut[1]=gcd(val1,n);
 vecOut[2]=gcd(val2,n);
 if( vecOut[1]==1 | vecOut[2] == 1, repeat=1;);
 a++; 

 );
 return(vecOut);

}


/*
**
**    log2(n)  
**
**    It works well with integer values n
**
*/
{log2(n) = local (k, i, nq, extrsup, stack);

  k=n;
  i=0;
  extrsup=0;
  stack=0;
  nq=0;
  
 /* if n isn't a power of 2 -> nq = 1 */
  
  if(!ispower(n,2), nq = 1;);
  
  
  while(k!=0,
    k=shift(k,-1);
    i++;
  );
   
  
  if(2^i > n, nq=1;);

  if( nq==1,
      extrsup = i;
      while( 2^extrsup > n, 
             extrsup=extrsup-0.00001; 
           );
      i=extrsup+1;   
   );
   stack=ceil(i-1);  /* This for precision: for example log2(32) = 5 */
   if( stack > i - 1 + 0.00001, stack=i-1;);
   return(stack);
}

/*
**    logb(n,b)  
*/
{logb(n,b) = local (k, l, m);
  
   k=log2(n);
   l=log2(b);
   m=1. * k/l;
   p=floor(m);
   if( m - p < 0.00001, m=p;); 
   return(m);
}

/*
** loge(n)  (Nepero)
*/
{loge(n) = local (k, extrsup, m);
  
   k=log2(n);

   extrsup = 2; /* 2 ^ 2 = 4 */
   while( 2^ extrsup > 2.71828182 , 
             extrsup=extrsup-0.00001; 
         );
   
   m = 1. * k/extrsup;

   return(m);
}

/*
**
**    isnCollatz(n)  
**
**    It returns 1 if n is a Collatz's number (2^k-1/3)
**               0 if no
*/
{isnCollatz(n) = local (value);

 value=0;

 appo=n*3+1;  /* It 's good the log2 beacuse k is always even */

 if( frac(log2(appo)) == 0.0, value = 1 ); 
 
 return(value);
}
/*
**
**    isnBizzarre(n)  
**
**    It returns 1 if n is a Bizzarre's number (2^k-4/4)
**               0 if no
*/
{isnBizzarre(n) = local (appo, value, appo2);

 value=0;
 appo=n*4+4;
 appo2=log2(appo);
 if(frac(appo2)==0.0, value=1);

 return(value);
}


/*
**    maxk2k(n)
**    It returns 2^k less than n or equal to n
**
*/

{maxk2k(n) = local (k);

 k=0;

 while( 2^k <= n, k++;);

 return(k-1);

}

/*
**    f2k(n)    
**    It returns the number as sum of power of 2
**
*/

{f2k(n) = local (k1,k2,k3,d1,d2,d3,s="");

 k1=maxk2k(n);
 d1=n-(2^k1);
 s=concat(s,"2^");
 s=concat(s,k1);
 
 while( d1>0,

        k2=maxk2k(d1);
        d2=d1-(2^k2);
        if( k2>=0, s=concat(s,"+2^"); s=concat(s,k2););

        k3=maxk2k(d2);
        d3=d2-(2^k3);

        d1=d3;

        if( k3>=0, s=concat(s,"+2^"); s=concat(s,k3););

        if(d3==1, s=concat(s,"+");s=concat(s,"2^0"); d1=0;);

 );
 return(s);

}

/*
**    is4np1(n)    
**    It returns if the number is a form 4n+1
**
*/
{is4np1(n) = local (v, t=0);

  v=(n-1)%4;
  if(v==0, t=1;);
  if(v!=0, t=0;);
  return(t);  
 
}


/*
**    isgc(n)    
**    It returns if the number is a form glide-complex
**
*/
{isgc(n) = local (s="", s2="", t=0, t1=0, t2=0, i=0, l=0);
  
 if( n%2 == 0, return(0);); /* ONLY ODD */

 appo=log2(n);

 if( frac(appo) == 0.0, return(0);); /* no power of 2 */
  
  s=f2k(n);
  l=length(s);
  s2=Vec(s);

  if( isnCollatz(n)==0 & isnBizzarre(n)==0,  

  for(i=1,l-1, 
      if(s2[i]=="^" & s2[i+1]=="4" & s2[i+2]=="+", t1=1; i=l;);
     );

  for(i=1,l-1, 
      if(s2[i]=="^" & s2[i+1]=="2" & s2[i+2]=="+", t2=1; i=l;);
     );
  );

  if(t1==0 & t2==0, t=1;);  /* If without 2^2 and 2^4 */

  if(t1==1 & t1==1, t=0;);  /* with 2^2 and 2^4 */

  if( isnCollatz(n)==1 | isnBizzarre(n)==1, t=0);  /* No Bizzarre's or Collatz's numbers */ 

  return(t);  

}

/*
**    is4np3(n)    
**    It returns if the number is a form 4n+3
**
*/
{is4np3(n) = local (v,t=0);
  
  v=(n-3)%4;
  if(v==0, t=1;);
  if(v!=0, t=0;);
  return(t);

}



/*
**    detectCnum(n)    
**    It returns the type of the number 
**
*/
{detectCnum(n) = local (s="", appo=0, appo1=0, appo2=0, appo3=0);
  
 s=concat(s,n);
 s=concat(s," = ");
 s = concat(s, f2k(n));
 s=concat(s," ");

 if( n%2==0, s=concat(s, "- even "););
 if( n%2!=0, s=concat(s, "- odd "););

 if( is4np3(n) == 1, s=concat(s,"- 4n+3 ");); 
 if( is4np1(n) == 1, s=concat(s,"- 4n+1 "););

 appo1=G1(n);
 appo2=G1(n+1);
 appo3=G1(n-1);


 /* if isoglide to a even */
 if( appo1==appo2, s=concat(s,"- G(");s=concat(s,n);s=concat(s,")=G(");s=concat(s,n+1);s=concat(s,") "););
 


 if( appo1==appo3, s=concat(s,"- G(");s=concat(s,n);s=concat(s,")=G(");s=concat(s,n-1);s=concat(s,") "););


 appo1=G1(n);
 appo2=G1(appo1);

 if(appo2==n,  s=concat(s," - glide-inverse ");); 
 
 if( isgc(n) == 1, s=concat(s," - glide-complex ");); 

 if( isnCollatz(n) == 1, s=concat(s,"- Collatz's number: ");
     appo = 3*n+1;
     appo2= log2(appo);
     s=concat(s,"(2^");
     s=concat(s,appo2);
     s=concat(s,"-1)/3 K+1=");
     s=concat(s,appo2);
     s=concat(s,"+1");


 );


 if( isnBizzarre(n) == 1, s=concat(s," - Bizzarre's number: ");
     appo = 4*n+4;
     appo2= log2(appo);
     s=concat(s,"(2^");
     s=concat(s,appo2);
     s=concat(s,"-4)/4");
 );

 appo=log2(n);  /* if power 2^ */
 if( frac(appo) == 0.0, s=concat(s, "- 2^"); s=concat(s, appo);); 
 

 appo=logb(n,3); /* if power 3^ */
 if( frac(appo) == 0.0, 
     s=concat(s, "- 3^"); 
     s=concat(s, appo);
     s=concat(s, ": G(");
     s=concat(s, n);
     s=concat(s, ")=G(");
     if(appo%2==0, appo2=(3*n+1)/4; s=concat(s,appo2); s=concat(s,")+3"); 
     );
     if(appo%2!=0, appo2=(3*n+1)/2; s=concat(s,appo2); s=concat(s,")+2");
     );
 ); 


 s=concat(s," - G(");
 s=concat(s,n);
 s=concat(s,")=");

 appo=0;
 appo=G1(n);
 s=concat(s,appo); 
 s=concat(s," "); 
 s=concat(s, isopm(n));
 
 G(n);

 return(s); 
}

/* 
*
*   gbiuniv(n1,n2)
*   It searches G(a)=b e G(b)=a
*
*/

{gbiuniv(n1,n2)=local(s="",i=1, appo1=0, appo2=0);

  for(i=n1,n2,

     appo1=G1(i);
     appo2=G1(appo1);

     if(appo2==i, 
        s=concat(s," G("); 
        s=concat(s,i); 
        s=concat(s,")=");
        s=concat(s,appo2);
        print(s);
        s=""

     );
     appo1=0;
     appo2=0;

  );
  
   
}

/* 
*
*   gplus(n, n1,n2)
*   It searches G(n1)=n1+n from n1 to n2
*
*/

{gplus(n, n1,n2,p=0)=local(s="",i=1, appo=0);

  for(i=n1,n2,

     appo=G1(i);

     if(appo==i+n, 
        s=concat(s," G("); 
        s=concat(s,i); 
        s=concat(s,")=");
        s=concat(s,appo);
        s=concat(s," ");
        if( p == 0, print(s); s="";);

     );
     appo=0;

  );
  if( p==1, return(s););  
   
}
/* 
*
*   gminus(n, n1,n2)
*   It searches G(n1)=n1-n from n1 to n2
*
*/

{gminus(n, n1,n2,p=0)=local(s="",i=1, appo=0);

  for(i=n1,n2,

     appo=G1(i);

     if(appo==i-n, 
        s=concat(s," G("); 
        s=concat(s,i); 
        s=concat(s,")=");
        s=concat(s,appo);
        s=concat(s," ");

        if( p == 0, print(s); s="";);

     );
     appo=0;

  );
  if( p==1, return(s););  
  
   
}


/* 
*
*   isog(n, n1,n2)
*   It searches G(n)=G(x) for x from n1 to n2
*
*/

{isog(n, n1,n2,p=0)=local(s="",i=1, appo1=0, appo2=0);

  appo1=G1(n);
  s=concat(s," G("); 
  s=concat(s,n); 
  s=concat(s,")");
  if( p == 0, print(s); s="";);

  for(i=n1,n2,

     appo2=G1(i);

     if(appo1==appo2 & i != n, 
        s=concat(s," = G(");
        s=concat(s,i);
        s=concat(s,") ");
 
        if( p == 0, print(s); s="";);


     );
     appo2=0;

  );
  if( p==1, return(s););  
   
}

/* 
*
*   isopm(n)
*   It searches the number isopath max
*
*/

{isopm(n)=local(s="",appo1=0);

  if( n%2 != 0,
      appo1=2*(3*n+1);
      s=concat(s," path max "); 
      s=concat(s," = G(");
      s=concat(s,appo1); 
      s=concat(s,") ");

  );
  return(s);  
   
}



/* 
*
*   getMatLevy(n)
*   It returns a Levy's Matrix
*   of numeric progression p+2q
*   where p and q are prime numbers
*   and n is nth prime number.
*
*/

{getMatLevy(n)=local(p=0,q=0,i=0,j=0);


     mat=matrix(n,n);
     
     for(p=2,n+1,

         for( q=2,n+1,
              i=p-1;
              j=q-1;
              mat[i,j] = primes(n+1)[p] + 2 * primes(n+1)[q];
         );  
      
     );
     return(mat);  
}



/* 
*
*   getCharMatLevy(n)
*   It returns The Charateristics Numbers (sums of the last two lines of Levy's matrix) and their difference 
*   
*
*/

{getCharMatLevy(n)=local(p=0,q=0,i=0,j=0);

     r1=0;
     r2=0;
     c1=0;
     c2=0;

     mat=matrix(n,n);
     
     for(p=2,n+1,

         for( q=2,n+1,
              i=p-1;
              j=q-1;
              mat[i,j] = primes(n+1)[p] + 2 * primes(n+1)[q];
         );  
      
     );

     if( n==2, m=1;);
     if( n>2,  m=n-1;);
     
     for( j=1,n,
          i=m;
          r1 = r1 + mat[i,j];
     );

     for( j=1,n,
          i=m+1;
          r2 = r2 + mat[i,j];
     );

    for(  i=1,n,
          j=m;
          c1 = c1 + mat[i,j];
     );

     for( i=1,n,
          j=m+1;
          c2 = c2 + mat[i,j];
     );

     print("r1 = ", r1);
     print("r2 = ", r2);
     print("c1 = ", c1);
     print("c2 = ", c2);
     print("r2-r1 = ", r2-r1);
     print("c2-c1 = ", c2-c1);
}
/* 
*
*   getMsngLevy(n)
*   It returns missing odd numbers in Levy's Matrix
*   n is nth prime number.
*
*/

{getMsngLevy(n)=local();


     p=0; q=0; exist=0; k=0; i=0; j=0; mx=0; 

     mat=matrix(n,n);
     c=n+1;

     for(p=2,c, for( q=2,c, i=p-1; j=q-1; mat[i,j] = primes(n+1)[p] + 2 * primes(n+1)[q]; );  );
  
  
     mx = primes(n+1)[n+1];
 
     m = (3*mx - 1)/2;

     d = 0;

     print(" ");
     print1( "Prime numbers are: " );

     for( i=2,n+1, print1(primes(n+1)[i]); print1(" ");  );  

     print1( "\nMissing Odd numbers: " );
     for(k=4,m,
         i=1;
         j=1;
         exist = 0;  
        
         d = 2 * k + 1;

         for(i=1,n,

              for( j=1,n, if(mat[i,j] ==  d, exist = 1;); );  
         );

         if( exist == 0, print1(d); print1(" "); );
      );

     getCharMatLevy(n);
}

/* 
*
*   getDetLevy(n)
*   It returns the determinant of a Levy's Matrix
*   of numeric progression p+2q
*   where p and q are prime numbers
*   and n is nth prime number.
*
*/

{getDetLevy(n)=local(p=0,q=0,i=0,j=0);


     mat=matrix(n,n);
     val=0;

     for(p=2,n+1,

         for( q=2,n+1,
              i=p-1;
              j=q-1;
              mat[i,j] = primes(n+1)[p] + 2 * primes(n+1)[q];
         );  
      
     );
     val = matdet(mat);
     return(val);  
}


/* 
*
*   listSigma(n)
*   It returns the list of sigma(n) (sum of divisors) from n1 to n2 
*
*/

{listSigma(n1,n2)=local(i=0);

     for(i=n1,n2,

        print(" ", sigma(i));      
     );
}



/* 
*
*   listIntRobin(n)
*   It returns the list of numbers which sigma(n)<exp(gamma) * n * ln(ln(n)) 
*   from n1 to n2
*
*/

{listIntRobin(n1,n2)=local(i=0);

     if( n1 == 0, n1=7;);
     c=1.781062;

     for(i=n1,n2,
        if( sigma(i) < c*i*log(log(i)), 
            print1(" ", i);      
        );
     );
}


/* 
*
*   listRatioRobin(n)
*   It returns the list of numbers with r=(exp(gamma) * n * ln(ln(n)))/sigma(n)
*   from n1 to n2
*
*/

{listRatioRobin(n1,n2)=local(i=0, r=0);

     if( n2==0, n2=n1;);
     c=1.781062;

     for(i=n1,n2,
            r= (c*i*log(log(i)))/sigma(i); 
            print(" ");
            print1("#", i);
            print1(" r= ", r);
     );
}