use"realroot";
//use"linalg";
//use"symbolix";
//include "mubase.mmx";
//R:=QQ['x];
//a:=R<<"x^5+x*x+1";
//b:=R<<"x^4+x*x+1";
//c:=R<<"x^3-x+1";
//d:=R<<"0";
POMR==> Polynomial MonomialSparseRing Rational;


degree(f,n:Int)==
{l:=#coefficients(f,n)-1;
 return l;
};

EEA(f:POMR,g:POMR): Vector Generic==
{
  zero:POMR := f-f;
  r0:POMR:=f;  r1:POMR:=g;// mmout<<r1<<"\n";
  s0:POMR:=zero+1;  s1:POMR:=zero;
  t0:POMR:=zero;  t1:POMR:=zero+1;
  if degree r1=0 and eval(r1,[0])=0 then return[r0,s0,t0];
  if degree r0=0 and eval(r0,[0])=0 then return[r1,s1,t1];
  while r1!=zero do {   
    q1:POMR:=r0/r1;
    r2:POMR:=r0-q1*r1;
    s2:POMR:=s0-q1*s1;
    t2:POMR:=t0-q1*t1;
    r0:=r1;
    r1:=r2;
    s0:=s1;
    s1:=s2;
    t0:=t1;
    t1:=t2;
  }
//}
 return [r0,s0,t0];
};

gcd(a:POMR,b:POMR):POMR=={
  f:POMR:=a;
  g:POMR:=b;
 // zero:POMR:=f-f;
  H:Vector Generic:=EEA(f,g);
  if degree a =0 and degree b=0 and eval(a,[0])=0 and eval(b,[0])=0 then return(a);
  return(H[0]/(coefficients H[0])[0]);
};

gcd(L: Vector Generic)=={
  n:Int:=#L;
  H:POMR:=gcd(L[0],L[1]);
  for i in 2..n do
    H:=gcd(H,L[i]);
  return(H);
};


GEEA(f,g,h)==
{Z:=  identitymatrix(3);
  s1:=f; s2:=g; s3:=h;
  zero:=f-f;
  while r!=zero do
  {M:=MatrixEuclide(f,g,h);
    q0:=s1/s2;
    r0:=rem(s1,s2);
    q1:=r0/s3;
    r:=s1-q0*s2-q1*s3;
    Z:=M*Z;
    s1:=s2;
    s2:=s3;
    s3:=r;
     mmout <<Z<<"\n";
    mmout <<s1<<","<<s2<<","<<s3<<"\n"; 
  }
  return(Z);
};