/********************************************************************

How to use this gadget:

Attach("s.m");

R<t> := FieldOfFractions(PolynomialRing(RationalField(),1));
 P := [<[t,0],2>, <[0,0],2>];

V, M, A := S(P,3);

********************************************************************/

intrinsic S(P::SeqEnum, d::RngIntElt) -> .
{Comment (without this, it's a syntax error).
P is a sequence of pairs <[x,y], m>. 
}
   require #P gt 0 : "Argument 1 must have cardinality at least 1.";
   require d ge 0 : "Argumenht 2 must be nonnegative.";
  
   K := Parent(P[1][1][1]);
   n := Binomial(d+2,2);
   R := PolynomialRing(K,n);
   AssignNames(~R,[Sprintf("a%o",i) : i in [1..n]]);
   S<X,Y> := PolynomialRing(R,2);
   mon := &cat[[m : m in MonomialsOfDegree(S,e)] : e in [0..d]];
   f := &+[R.i*mon[i] : i in [1..n]];
   nrels := &+[Binomial(p[2]+1,2) : p in P];
   A := RMatrixSpace(K,n,nrels)!0;
   c := 1;   // current column.
   for i in [1..#P] do
      Q := P[i][1];
      m := P[i][2];
      v := { [] };
      for deg in [0..m-1] do
         /* Compute the linear conditions coming from partials of degree deg
            and insert them in the matrix A as columns. */
          for D in v do
             g := Derivative(f,D);
             // L is a linear form that expresses the condition 
             // that D(f) vanishes at the point Q.
             L := Evaluate(g,Q);
             for r in [1..Nrows(A)] do 
                A[r,c] := MonomialCoefficient(L,R.r);
             end for;
             c := c + 1;
          end for;
          v := { Sort(D cat [1]) : D in v } join { Sort(D cat [2]) : D in v } ;
      end for;
   end for;
   V := Kernel(A);
   M := [&+[mon[i]*v[i] : i in [1..#mon]] : v in Basis(V)];
   return V, M, A;
end intrinsic;


