//
proc nFactorial(int n)
{
 // Compute n factorial for a given integer n
    if(n==0)
    {
        return(1);
    }
    else
    {
        bigint k=1;
        int i=1;
        while(i<=n)
        {
            k=k*i;
            //k;
            i++;
        }
        return(k);
     }
}

proc nFact(int n)
{
   if(n<0)
   {
       // "------ WARNING ----------"
     ERROR(" Warning: n must be greater or equal to 0 ");
   }
   bigint k =1;
   while(n >=2)
   {
      k = k*n;
      n = n-1;
   }
   return(k);
}

proc Trace(intmat m)
{
    int k;
    for(int j=1; j<=nrows(m); j++)
    {
        k= k+m[j,j];
    }
       return(k);
}


proc Gcd(int n, int m)
{
    int a, b;
    a = m; 
    b = n;
    if(a == 0 and  b == 0)
    {
        return(0);
    }
    else
    {
        int r;
        while( b != 0)
        {
            r=a%b;
            a=b;
            b=r;
        }
        if(a<0)
        {
            return(-a);
        }
        return(a);
    }
}

// compute a list of primes less than or equal to an integer n
proc listprime(int n)
{
    if(n==2)
    {
        return(list(2));
    }
    if(n<=1)
    {
         return(list());// list() is empty list
    }
    else
    {
         // here n> 2
         list T=2..n;
         intvec v=T[1];
         list l=v[1..nrows(v)]; // l = 2,3, 4,5,6, ...
         list m;
         int i=1;
         int j; // j = 0
         while(1)
         {
              if(l[i]^2>n)
              {
                  return(l);
              }
              for(j=i+1;j<=size(l);j++)
              {
                   if(l[j]%l[i]!=0)
                   {
                        m= m+list(l[j]);
                   }
              }
              m=list(l[1..i])+m;
              l=m;
              m=list(); // empty list 
              i=i+1;
         }
     }
}                          

proc list_prime(int n)
{
    int k=n;
    list l;
    while(k>1)
    {
        k=prime(k)-1;
        l=l+list(k+1);
    }
    return(l);
}

proc lprime(int n)
{
    list l;
    for(int i = 2; i<=n; i++)
    {
        if(prime(i) == i)
        {
            l = l + list(i);
        }
    }
    return(l);
}

proc extGCD(bigint a,  bigint b)
{
    bigint r0,r1,s0,s1,t0,t1,q,r2,s2,t2;
    r0=a; r1=b; s0=1; t0=0;
    s1=t0; t1=s0;
    //list l;
    while(r1!=0)
    {
        //l=division(r0,r1);
        //r2 = r0%r1;
        q = r0 div r1;
        r2 = r0 - r1*q;
        //q= l[1][1,1];
        //r2= l[2][1];
        s2=s0-q*s1;
        t2=t0-q*t1;
        r0=r1; r1=r2;
        s0=s1; s1=s2;
        t0=t1; t1=t2;
        //"=========================";
        //(r0,s0,t0);
    }
    //return(r0,s0,t0);
    return("Dereje")
}                         
proc extGCDpoly(poly a,  poly b)
{
    poly r0,r1,s0,s1,t0,t1,q,r2,s2,t2;
    r0=a; r1=b; s0=1; t0=0;
    s1=t0; t1=s0;
    list l;
    while(r1!=0)
    {
        l = division(r0,r1);
        q = l[1][1,1];
        r2 = l[2][1];
        s2 = s0-q*s1;
        t2 = t0-q*t1;
        r0 = r1; r1 = r2;
        s0 = s1; s1 = s2;
        t0 = t1; t1 = t2;
        "=========================";
        (r0,s0,t0);
    }
    return(r0, s0, t0);
}          

proc normN(poly f)
{
    if(f!=0)
    {
        return(f/leadcoef(f));
    }
    else
    {
        return(lu(f));
    }
}                


proc lu(poly f)
{
    if(f!=0)
    {
        return(leadcoef(f));
    }
    else
    {
          return(poly(1));
    }
}           


proc GcdpolyN(poly a, poly b)
{
    poly h0,h1,h2,r0,r1,s0,s1,t0,t1,q,r2,s2,t2;
    h0=lu(a); h1=lu(b);
    r0=normN(a); r1=normN(b); s0=1/h0;t0=0;
    s1=t0; t1=1/h1;
    list l;
    while(r1!=0)
    {
        //"===========";
        //(r0,s0,t0);
        l=division(r0,r1);// divide ro by r1
        q=l[1][1,1];// q = ro quo r1
        r2=l[2][1];// r2 = ro rem r1
        h2=lu(r2);
        s2=(s0-q*s1)/h2;
        t2=(t0-q*t1)/h2;
        r0=r1; r1=r2;
        s0=s1; s1=s2;
        t0=t1; t1=t2;
    }
    return(cleardenom(r0)); //list (r0,s0,t0));
}