include"mubase.mmx";
include"pencilregular.mmx";
//R1:=QQ['x,'y,'z,'w];
//var1:=R1<<"x";
//var2:=R1<<"y";
//var3:=R1<<"z";
//var4:=R1<<"w";
//varsur:=[var1,var2,var3];
//varcur:=[var1,var2];
//varimplcur:=[var1,var2,var3];
//varimplsur:=[var1,var2,var3,var4];
//f0:=R1<<"x*z^2";
//f1:=R1<<"y^2*x+y^2*z";
//f2:=R1<<"x^2*y+x*y*z";
//f3:=R1<<"x*y*z+y*z^2";
//g0:=R1<<"x^2+y^2+z^2";
//g1:=R1<<"2*x*z";
//g2:=R1<<"2*x*y";
//g3:=R1<<"x^2-y^2-z^2";
//cur0:=R1<<"y^3";
//cur1:=R1<<"y^2*x";
//cur2:=R1<<"y*x^2";
//cur3:=R1<<"x^3";
//p0:= R1<<"x^2";
//p1:= R1<<"x*y";
//p2:= R1<<"y^2";
//q0:=R1<<"x^3+y^3";
//q1:=R1<<"x^2*y";
//q2:=R1<<"x^3+x*y^2+y^3";

//a11:=R1<<"3*x^2*y";
//a21:=R1<<"-3*x^3+3*x*y^2";
//a31:=R1<<"x*y^2+y^3";
//a41:=R1<<"y^3";
//stein1:=R1<<"x^2+y^2+z^2";
//stein2:=R1<<"y*z";
//stein3:=R1<<"x*z";
//stein4:=R1<<"x*y";
//stein:=[stein1,stein2,stein3,stein4];
//sphere:=[g0,g1,g2,g3];
//curstein:= [R1<<"x^3+x*y^2+y^3",cur2,cur1,cur0];
matrix_multiply_curve(p: POMR,var: Vector Generic,n:Int,d:Int)==
{ f:POMR:=(var[0]+var[1])^n; 
   A:Vector Generic:= monomials f;
   a1:Int:=#A;
   g:POMR:=(var[0]+var[1])^(n+d);
   B:Vector Generic:= monomials g;
   b1:Int:=#B;
   Z:Matrix Rational:=zero_matrix_rational(b1,a1);
     
   for i in 0..a1 do
     {
   C:POMR:=A[i]*p;
   D:Vector Generic:= monomials C;
   E:Vector Generic:= coefficients C;
   for j in 0..#E do
     for k in 0..b1 do
       if D[j]=B[k] then Z[k,i]:=E[j];
   }
  return(Z); 
};
////f:=[f0: POMR,f1: POMR,f2: POMR,f3: POMR], var:=[var1: POMR,var2: POMR,var3: POMR] of f

matrix_koszul_curve_plane(f: Vector Generic,var: Vector Generic)==
 {    n:Int:= degree f[0]+degree f[1]-1;
      L:Matrix Rational:= matrix_multiply_curve(f[0],var,n-degree f[0],degree f[0]);
      M:Matrix Rational:= matrix_multiply_curve(f[1],var,n-degree f[1],degree f[1]);
      N:Matrix Rational:= matrix_multiply_curve(f[2],var,n-degree f[2],degree f[2]);
      Z:Matrix Rational:=horizontal_join(horizontal_join(L,M),N);
   return(Z);
 };
///f:=[f0,f1,f2,f3], var1:=[var11,var12,var13], variable of f; var2:=[var20,var21,var22,var23];

//matrix_rep_curve_plane(f:Vector Generic,varcur: Vector Generic,varimpl: Vector Generic)== {
 //  A:Matrix Rational:=matrix_koszul_curve_plane(f,varcur);
  // C:Matrix Rational:= transpose ker A; 
  // d:Int:=degree f[0];
  //  m:Int:=(degree f[0])*(degree f[0]+1)div 2:>Int;
  // Z:=zero_matrix_generic(d,rows C);
  // for i in 0..d do
    // for j in 0..rows C do
     //  Z[i,j]:=varimpl[0]*C[j,i]+varimpl[1]*C[j,i+d]+varimpl[2]*C[j,i+2*d];
  // return(Z);
 //};
 

 matrix_multiply_surface(p: POMR,var: Vector Generic,n:Int,d:Int)==
{  f:POMR:=(var[0]+var[1]+var[2])^n; 
   A:Vector Generic:= monomials f;
   a1:Int:=#A;
   g:POMR:=(var[0]+var[1]+var[2])^(n+d);
   B:Vector Generic:= monomials g;
   b1:Int:=#B;
   Z:Matrix Rational:=zero_matrix_rational(b1,a1);
     
   for i in 0..a1 do
     {
   C:POMR:=A[i]*p;
   D:Vector Generic:= monomials C;
   E:Vector Generic:= coefficients C;
   for j in 0..#E do
     for k in 0..b1 do
       if D[j]=B[k] then Z[k,i]:=E[j];
   }
  return(Z); 
};
////f:=[f0: POMR,f1: POMR,f2: POMR,f3: POMR], var:=[var1: POMR,var2: POMR,var3: POMR] of f
 matrix_koszul_surface(f: Vector Generic,var: Vector Generic)==
 {  n:Int:= degree f[0]+degree f[1]+degree f[2]-2;
    L:Matrix Rational:= matrix_multiply_surface(f[0],var,n-degree f[0],degree f[0]);
    M:Matrix Rational:= matrix_multiply_surface(f[1],var,n-degree f[1],degree f[1]);
    N:Matrix Rational:= matrix_multiply_surface(f[2],var,n-degree f[2],degree f[2]);
    K:Matrix Rational:= matrix_multiply_surface(f[3],var,n-degree f[3],degree f[3]);
    Z:Matrix Rational:=horizontal_join(horizontal_join(horizontal_join(L,M),N),K);
   return(Z);
 };
///f:=[f0,f1,f2,f3], var1:=[var11,var12,var13], variable of f; var2:=[var20,var21,var22,var23];

 matrix_rep_surface(f:Vector Generic,varsur: Vector Generic,varimpl: Vector Generic)== {
   A:Matrix Rational:=matrix_koszul_surface(f,varsur);
   // mmout<< A<<"\n";
  // B:Matrix Rational:=ker A;
   C:Matrix Rational:= transpose ker A; 
  // n:Int:=rows C;
   // d:Int:=degree f[0];
   //mmout<<rows C<<"\n";
   m:Int:=degree f[0]*(2*degree f[0]-1);
     //div 2:>Int;
   Z:Matrix Generic:=zero_matrix_generic(m,rows C);
   for i in 0..m do
     for j in 0..rows C do
       Z[i,j]:=varimpl[0]*C[j,i]+varimpl[1]*C[j,i+m]+varimpl[2]*C[j,i+2*m]+varimpl[3]*C[j,i+3*m];
   return(Z);
 };

matrix_vector_rep(M: Matrix Generic, cur: Vector Generic,n:Int): Matrix Generic=={
   Z: Matrix Generic:=[as_polynomial_POMR(0,cur[0])|i in 0..n||j in 0..1];  
   for i in 0..n do
     for j in 0..columns M do
       Z[0,i]:=Z[0,i]+cur[j]*M[0,j][i];
   return(Z);
};

matrix_rep_curve_plane_homogeneous(cur: Vector Generic,varcur: Vector Generic,curimpl: Vector Generic)==
{curve:Vector Generic:=[];
   for i in 0..#cur do
    { a:POMR:=as_polynomial_POMR(eval(cur[i],[varcur[0],1]),varcur[0]);
     curve:=cons(a,curve)};
   curve:=reverse curve;
   g:Vector Generic:=zero_vector_generic(#cur);
   for i in 0..#cur do
     g[i]:=as_polynomial_rational curve[i];
 L:Vector Generic:=mubase g;
    A: Matrix Generic:=transpose(matrix_vector_rep(L[0],curimpl,degree L[0]+ degree L[1]));
  for i in 1..degree L[1] do
    // A:=horizontal_join(A,transpose matrix_vector_rep(polynomial(0,1)^i*L[0],curimpl,degree L[0]+degree L[1]));
     A:=horizontal_join(A,transpose matrix_vector_rep(matrixp(L[0],polynomial(0,1)^i),curimpl,degree L[0]+degree L[1]));
  //mmout<<A<<"\n";
  B: Matrix Generic:=transpose matrix_vector_rep(L[1],curimpl,degree L[0]+degree L[1]);
  for i in 1..degree L[0] do
    // B:=horizontal_join(B,transpose matrix_vector_rep(polynomial(0,1)^i*L[1],curimpl,degree L[0]+degree L[1]));
    B:=horizontal_join(B,transpose matrix_vector_rep(matrixp(L[1],polynomial(0,1)^i),curimpl,degree L[0]+degree L[1]));
    C: Matrix Generic:=horizontal_join(A,B);
  return(C);
};

matrix_rep_curve_plane(cur: Vector Generic,curimpl: Vector Generic)==
{ //curve:Vector Generic:=[];
   //for i in 0..#cur do
   // { a:POMR:=as_polynomial_POMR(eval(cur[i],[varcur[0],1]),varcur[0]);
    // curve:=cons(a,curve)};
   //curve:=reverse curve;
   g:Vector Generic:=zero_vector_generic(#cur);
   for i in 0..#cur do
     g[i]:=as_polynomial_rational cur[i];
 L:Vector Generic:=mubase g;
    A: Matrix Generic:=transpose(matrix_vector_rep(L[0],curimpl,degree L[0]+ degree L[1]));
  for i in 1..degree L[1] do
    // A:=horizontal_join(A,transpose matrix_vector_rep(polynomial(0,1)^i*L[0],curimpl,degree L[0]+degree L[1]));
     A:=horizontal_join(A,transpose matrix_vector_rep(matrixp(L[0],polynomial(0,1)^i),curimpl,degree L[0]+degree L[1]));
  //mmout<<A<<"\n";
  B: Matrix Generic:=transpose matrix_vector_rep(L[1],curimpl,degree L[0]+degree L[1]);
  for i in 1..degree L[0] do
    // B:=horizontal_join(B,transpose matrix_vector_rep(polynomial(0,1)^i*L[1],curimpl,degree L[0]+degree L[1]));
    B:=horizontal_join(B,transpose matrix_vector_rep(matrixp(L[1],polynomial(0,1)^i),curimpl,degree L[0]+degree L[1]));
    C: Matrix Generic:=horizontal_join(A,B);
  return(C);
};



matrix_rep_curve_space_homogeneous(cur: Vector Generic,varcur: Vector Generic,curimpl: Vector Generic)==
{curve:Vector Generic:=[];
   for i in 0..#cur do
    { a:POMR:=as_polynomial_POMR(eval(cur[i],[varcur[0],1]),varcur[0]);
     curve:=cons(a,curve)};
   curve:=reverse curve;
   g:Vector Generic:=zero_vector_generic(#cur);
   for i in 0..#cur do
     g[i]:=as_polynomial_rational curve[i];
 L:Vector Generic:=mubase g;
    A: Matrix Generic:=transpose(matrix_vector_rep(L[0],curimpl,degree L[0]+ degree L[1]));
  for i in 1..degree L[1] do
    //   A:=horizontal_join(A,transpose matrix_vector_rep(polynomial(0,1)^i*L[0],curimpl,degree L[0]+degree L[1]));
    A:=horizontal_join(A,transpose matrix_vector_rep(matrixp(L[0],polynomial(0,1)^i),curimpl,degree L[0]+degree L[1]));
  //mmout<<A<<"\n";
  B: Matrix Generic:=transpose matrix_vector_rep(L[1],curimpl,degree L[0]+degree L[1]);
  for i in 1..degree L[0] do
    // B:=horizontal_join(B,transpose matrix_vector_rep(polynomial(0,1)^i*L[1],curimpl,degree L[0]+degree L[1]));
    B:=horizontal_join(B,transpose matrix_vector_rep(matrixp(L[1],polynomial(0,1)^i),curimpl,degree L[0]+degree L[1]));
   // mmout<<B<<"\n";
   C: Matrix Generic:=transpose matrix_vector_rep(L[2],curimpl,degree L[0]+degree L[1]);
 
  for i in 1..(degree L[0]+degree L[1]-degree L[2]) do
    //  C:=horizontal_join(C,transpose matrix_vector_rep(polynomial(0,1)^i*L[2],curimpl,degree L[0]+degree L[1]));
    C:=horizontal_join(C,transpose matrix_vector_rep(matrixp(L[2],polynomial(0,1)^i),curimpl,degree L[0]+degree L[1]));
   // mmout<<C<<"\n";
  D: Matrix Generic:=horizontal_join(horizontal_join(A,B),C);
  return(D);
};



matrix_vector(M: Matrix Generic, cur: Vector Generic,n:Int): Matrix Generic=={
   Z: Matrix Generic:=zero_matrix_generic(1, n);
   for i in 0..n do
     for j in 0..columns M do
       Z[0,i]:=Z[0,i]+cur[j]*M[0,j][i];
   return(Z);
};


matrix_rep_curve_space(cur1: Vector Generic,curimpl: Vector Generic)=={
  L:Vector Generic:=mubase cur1; //mmout<<L<<"\n";
  A: Matrix Generic:=transpose(matrix_vector_rep(L[0],curimpl,degree L[0]+ degree L[1]));
  for i in 1..degree L[1] do
    //  A:=horizontal_join(A,transpose matrix_vector_rep(polynomial(0,1)^i*L[0],curimpl,degree L[0]+degree L[1]));
     A:=horizontal_join(A,transpose matrix_vector_rep(matrixp(L[0],polynomial(0,1)^i),curimpl,degree L[0]+degree L[1]));
  //mmout<<A<<"\n";
  B: Matrix Generic:=transpose matrix_vector_rep(L[1],curimpl,degree L[0]+degree L[1]);
  for i in 1..degree L[0] do
    // B:=horizontal_join(B,transpose matrix_vector_rep(polynomial(0,1)^i*L[1],curimpl,degree L[0]+degree L[1]));
    B:=horizontal_join(B,transpose matrix_vector_rep(matrixp(L[1],polynomial(0,1)^i),curimpl,degree L[0]+degree L[1]));
   // mmout<<B<<"\n";
   C: Matrix Generic:=transpose matrix_vector_rep(L[2],curimpl,degree L[0]+degree L[1]);
  for i in 1..(degree L[0]+degree L[1]-degree L[2]) do
    // C:=horizontal_join(C,transpose matrix_vector_rep(polynomial(0,1)^i*L[2],curimpl,degree L[0]+degree L[1]));
     C:=horizontal_join(C,transpose matrix_vector_rep(matrixp(L[2],polynomial(0,1)^i),curimpl,degree L[0]+degree L[1]));
   // mmout<<C<<"\n";
  D: Matrix Generic:=horizontal_join(horizontal_join(A,B),C);
  return(D);
};

 //// Use mu-basis

matrix_rep_intersection_curve_plane(cur1: Vector Generic,cur2: Vector Generic)=={
  L:Vector Generic:=mubase cur1; // mmout<<L<<"\n";
  A: Matrix Generic:=transpose(matrix_vector(L[0],cur2,degree L[0]+ degree L[1]));
  for i in 1..degree L[1] do
    // A:=horizontal_join(A,transpose matrix_vector(polynomial(0,1)^i*L[0],cur2,degree L[0]+degree L[1]));
    A:=horizontal_join(A,transpose matrix_vector(matrixp(L[0],polynomial(0,1)^i),cur2,degree L[0]+degree L[1]));
  B: Matrix Generic:=transpose matrix_vector(L[1],cur2,degree L[0]+degree L[1]);
  for i in 1..degree L[0] do
    //  B:=horizontal_join(B,transpose matrix_vector(polynomial(0,1)^i*L[1],cur2,degree L[0]+degree L[1]));
     B:=horizontal_join(B,transpose matrix_vector(matrixp(L[1],polynomial(0,1)^i),cur2,degree L[0]+degree L[1]));
  C: Matrix Generic:=horizontal_join(A,B);
  return(C);
};

intersection_curve_plane(cur1: Vector Generic, cur2: Vector Generic)=={
   curve1:Vector Generic:=zero_vector_generic(#cur1);
   for i in 0..#cur1 do
     curve1[i]:=as_polynomial_rational cur1[i];
   curve2:Vector Generic:=zero_vector_generic(#cur2);
   for i in 0..#cur2 do
     curve2[i]:=as_polynomial_rational cur2[i];
  M:Matrix Generic:= matrix_rep_intersection_curve_plane(curve1,curve2);
  L:Vector Generic:=generalized_eigenvalues M;
  ListPoint:Vector Generic:=[];
  for i in 0..#L do
    ListPoint:=cons([eval(curve2[0],L[i]),eval(curve2[1],L[i]),eval(curve2[2],L[i])],ListPoint);
  return(delete_dif ListPoint);
};

/// Type: POMR

matrix_rep_intersection_curve_plane_homogeneous(cur1: Vector Generic,varcur1: Vector Generic,cur2: Vector Generic,varcur2: Vector Generic)=={
  curve1:Vector Generic:=[];
   for i in 0..#cur1 do
    { a:POMR:=as_polynomial_POMR(eval(cur1[i],[varcur1[0],1]),varcur1[0]);
     curve1:=cons(a,curve1)};
   curve1:=reverse curve1;
   g:Vector Generic:=zero_vector_generic(#cur1);
   for i in 0..#cur1 do
     g[i]:=as_polynomial_rational curve1[i];
   L:Vector Generic:=mubase g;
   //mmout<<L<<"\n";
   curve2:Vector Generic:=[];
   for i in 0..#cur2 do
    { a:POMR:=as_polynomial_POMR(eval(cur2[i],[varcur2[0],1]),varcur2[0]);
     curve2:=cons(a,curve2)};
   curve2:=reverse curve2;
    h:Vector Generic:=zero_vector_generic(#cur2);
   for i in 0..#cur2 do
     h[i]:=as_polynomial_rational curve2[i];
  // mmout<<h<<"\n";
   A: Matrix Generic:=transpose(matrix_vector(L[0],h,degree L[0]+ degree L[1]));
  // mmout << A<<"\n";
  for i in 1..degree L[1] do
    //  A:=horizontal_join(A,transpose matrix_vector(polynomial(0,1)^i*L[0],h,degree L[0]+degree L[1]));
     A:=horizontal_join(A,transpose matrix_vector(matrixp(L[0],polynomial(0,1)^i),h,degree L[0]+degree L[1]));
 // mmout<<A<<"\n";
  B: Matrix Generic:=transpose matrix_vector(L[1],h,degree L[0]+degree L[1]);
 // mmout <<B<<"\n";
  for i in 1..degree L[0] do
    //B:=horizontal_join(B,transpose matrix_vector(polynomial(0,1)^i*L[1],h,degree L[0]+degree L[1]));
    B:=horizontal_join(B,transpose matrix_vector(matrixp(L[1],polynomial(0,1)^i),h,degree L[0]+degree L[1]));

  // mmout<<B<<"\n";
    M: Matrix Generic:=horizontal_join(A,B);  
  return(M); 
 };

intersection_curve_plane_homogeneous(cur1: Vector Generic,varcur1: Vector Generic,cur2: Vector Generic, varcur2:Vector Generic)=={

  M: Matrix Generic:=matrix_rep_intersection_curve_plane_homogeneous(cur1,varcur1,cur2,varcur2);
  L:Vector Generic:=generalized_eigenvalues M;
  ListPoint:Vector Generic:=[];
  
  curve2:Vector Generic:=[];
   for i in 0..#cur2 do
    { a:POMR:=as_polynomial_POMR(eval(cur2[i],[varcur2[0],1]),varcur2[0]);
     curve2:=cons(a,curve2)};
   curve2:=reverse curve2;
    h:Vector Generic:=zero_vector_generic(#cur2);
   for i in 0..#cur2 do
     h[i]:=as_polynomial_rational curve2[i];
   for i in 0..#L do
    ListPoint:=cons([eval(h[0],L[i]),eval(h[1],L[i]),eval(h[2],L[i])],ListPoint);
 return(delete_dif ListPoint);  
};

matrix_rep_intersection_curve_space(cur1: Vector Generic,cur2: Vector Generic)=={
  L:Vector Generic:=mubase cur1; //mmout<<L<<"\n";
  A: Matrix Generic:=transpose(matrix_vector(L[0],cur2,degree L[0]+ degree L[1]));
  for i in 1..degree L[1] do
    //  A:=horizontal_join(A,transpose matrix_vector(polynomial(0,1)^i*L[0],cur2,degree L[0]+degree L[1]));
     A:=horizontal_join(A,transpose matrix_vector(matrixp(L[0],polynomial(0,1)^i),cur2,degree L[0]+degree L[1]));
  //mmout<<A<<"\n";
  B: Matrix Generic:=transpose matrix_vector(L[1],cur2,degree L[0]+degree L[1]);
  for i in 1..degree L[0] do
    //   B:=horizontal_join(B,transpose matrix_vector(polynomial(0,1)^i*L[1],cur2,degree L[0]+degree L[1]));
    B:=horizontal_join(B,transpose matrix_vector(matrixp(L[1],polynomial(0,1)^i),cur2,degree L[0]+degree L[1]));
   // mmout<<B<<"\n";
   C: Matrix Generic:=transpose matrix_vector(L[2],cur2,degree L[0]+degree L[1]);
  for i in 1..(degree L[0]+degree L[1]-degree L[2]) do
    // C:=horizontal_join(C,transpose matrix_vector(polynomial(0,1)^i*L[2],cur2,degree L[0]+degree L[1]));
    C:=horizontal_join(C,transpose matrix_vector(matrixp(L[2],polynomial(0,1)^i),cur2,degree L[0]+degree L[1]));
   // mmout<<C<<"\n";
  D: Matrix Generic:=horizontal_join(horizontal_join(A,B),C);
  return(D);
};

intersection_curve_space(cur1: Vector Generic, cur2: Vector Generic)=={
  curve1:Vector Generic:=zero_vector_generic(#cur1);
   for i in 0..#cur1 do
     curve1[i]:=as_polynomial_rational cur1[i];
   curve2:Vector Generic:=zero_vector_generic(#cur2);
   for i in 0..#cur2 do
     curve2[i]:=as_polynomial_rational cur2[i];
  M:Matrix Generic:= matrix_rep_intersection_curve_space(curve1,curve2);
  L:Vector Generic:=generalized_eigenvalues M;
  ListPoint:Vector Generic:=[];
  for i in 0..#L do
    ListPoint:=cons([eval(curve2[0],L[i]),eval(curve2[1],L[i]),eval(curve2[2],L[i]),eval(curve2[3],L[i])],ListPoint);
 return(delete_dif ListPoint);  
};


intersection_curve_space_homogeneous(cur1: Vector Generic,varcur1: Vector Generic,cur2: Vector Generic, varcur2:Vector Generic)==
{curve1:Vector Generic:=[];
   for i in 0..#cur1 do
    { a:POMR:=as_polynomial_POMR(eval(cur1[i],[varcur1[0],1]),varcur1[0]);
     curve1:=cons(a,curve1)};
   curve1:=reverse curve1;
   g:Vector Generic:=zero_vector_generic(#cur1);
   for i in 0..#cur1 do
     g[i]:=as_polynomial_rational curve1[i];
   L:Vector Generic:=mubase g;
   
   curve2:Vector Generic:=[];
   for i in 0..#cur2 do
    { a:POMR:=as_polynomial_POMR(eval(cur2[i],[varcur2[0],1]),varcur2[0]);
     curve2:=cons(a,curve2)};
   curve2:=reverse curve2;
   h:Vector Generic:=zero_vector_generic(#cur2);
   for i in 0..#cur2 do
     h[i]:=as_polynomial_rational curve2[i];
    A: Matrix Generic:=transpose(matrix_vector(L[0],h,degree L[0]+ degree L[1]));
  for i in 1..degree L[1] do
   // A:=horizontal_join(A,transpose matrix_vector(polynomial(0,1)^i*L[0],h,degree L[0]+degree L[1]));
  A:=horizontal_join(A,transpose matrix_vector(matrixp(L[0],polynomial(0,1)^i),h,degree L[0]+degree L[1]));
  //mmout<<A<<"\n";
  B: Matrix Generic:=transpose matrix_vector(L[1],h,degree L[0]+degree L[1]);
  for i in 1..degree L[0] do
   // B:=horizontal_join(B,transpose matrix_vector(polynomial(0,1)^i*L[1],h,degree L[0]+degree L[1]));
   B:=horizontal_join(B,transpose matrix_vector(matrixp(L[1],polynomial(0,1)^i),h,degree L[0]+degree L[1]));
   // mmout<<B<<"\n";
   C: Matrix Generic:=transpose matrix_vector(L[2],h,degree L[0]+degree L[1]);
  for i in 1..(degree L[0]+degree L[1]-degree L[2]) do
    //  C:=horizontal_join(C,transpose matrix_vector(polynomial(0,1)^i*L[2],h,degree L[0]+degree L[1]));
     C:=horizontal_join(C,transpose matrix_vector(matrixp(L[2],polynomial(0,1)^i),h,degree L[0]+degree L[1]));
   // mmout<<C<<"\n";
  M: Matrix Generic:=horizontal_join(horizontal_join(A,B),C);
   L:Vector Generic:=generalized_eigenvalues M;
  ListPoint:Vector Generic:=[];
  for i in 0..#L do
    ListPoint:=cons([eval(h[0],L[i]),eval(h[1],L[i]),eval(h[2],L[i]),eval(h[3],L[i])],ListPoint);
 return(delete_dif ListPoint);  
};


//// var:=[var1,var2,var3] list of variables of polynomial; type var1,..: POMR

////f:=[f0: POMR,f1: POMR,f2: POMR,f3: POMR], var:=[var1: POMR,var2: POMR,var3: POMR] of f
 
 //f:=[f0: POMR, f1: POMR,f2: POMR,f3:POMR],var:=[var1: POMR,var2: POMR,var3: POMR], cur:=[cur0,cur1,cur2,cur3]
 matrix_rep_intersection_curve_surface(f: Vector Generic,varsur: Vector Generic,cur: Vector Generic): Matrix Generic=={
   A:Matrix Rational:=matrix_koszul_surface(f,varsur);
 //  B:Matrix Rational:=ker A;
   C:Matrix Rational:= transpose ker A; 
  // n:Int:=rows C;
  // d:Int:=degree f[0];
   m:Int:=(degree f[0])*(2*degree f[0]-1);
   //div 2:>Int;
   Z:Matrix Generic:=zero_matrix_generic(m,rows C);
   for i in 0..m do
     for j in 0..rows C do
       Z[i,j]:=as_polynomial_rational(cur[0])*C[j,i]+as_polynomial_rational(cur[1])*C[j,i+m]+as_polynomial_rational(cur[2])*C[j,i+2*m]+as_polynomial_rational(cur[3])*C[j,i+3*m];
   return(Z);
 };
 
intersection_curve_surface(sur: Vector Generic,varsur: Vector Generic, cur: Vector Generic)=={
 M:Matrix Generic:= matrix_rep_intersection_curve_surface(sur,varsur,cur);
 L:Vector Generic:=generalized_eigenvalues M;
 ListPoint:Vector Generic:=[];
  for i in 0..#L do
    ListPoint:=cons([eval(cur[0],L[i]),eval(cur[1],L[i]),eval(cur[2],L[i]),eval(cur[3],L[i])],ListPoint);
return(delete_dif ListPoint);
};

intersection_curve_surface_homogeneous(sur: Vector Generic,varsur: Vector Generic, cur: Vector Generic,varcur: Vector Generic)=={
  curve: Vector Generic:=[eval(cur[0],[varcur[0],1]),eval(cur[1],[varcur[0],1]),eval(cur[2],[varcur[0],1]),eval(cur[3],[varcur[0],1])]; 
 M:Matrix Generic:= matrix_rep_intersection_curve_surface(sur,varsur,curve); 
 L:Vector Generic:=generalized_eigenvalues M;
// mmout<<L<<"\n";
  curve2: Vector Generic:=[eval(cur[0],[1,0]),eval(cur[1],[1,0]),eval(cur[2],[1,0]),eval(cur[3],[1,0])]; 
 M2:Matrix Generic:= matrix_rep_intersection_curve_surface(sur,varsur,curve2); 

 if rank M2 < min(rows M2, columns M2) then ListPoint:Vector Generic:=[curve2];
 else ListPoint:Vector Generic:=[];
 //mmout<< ListPoint<<"\n";
  for i in 0..#L do
    ListPoint:=cons([eval(as_polynomial_rational(curve[0]),L[i]),eval(as_polynomial_rational(curve[1]),L[i]),eval(as_polynomial_rational(curve[2]),L[i]),eval(as_polynomial_rational(curve[3]),L[i])],ListPoint);
return(delete_dif ListPoint);
};