MatsuoAlgebra:=function(R,M,c,h)
 D:=[1..Nrows(M)];
 nl:=Ncols(M);
 Mt:=Transpose(M);
 MMt:=M*Mt;
 L:=[ Support(Mt[k]) : k in [1..nl] ];
 V:=VectorSpace(Rationals(),#D);
 cv:=func< s | CharacteristicVector(V,s) >;
 line:=func< i,j | [ k : k in [1..nl] | {i,j} subset L[k] ] >;
 thirdpt:=func< i,j | L[line(i,j)[1]] diff {i,j} >;
 // structure constants
 scv:=function(i,j)
  if i eq j then return 2*cv({i});
  else L:=line(i,j);
   if L eq [] then return cv({});
   else return 1/2*h*(cv({i,j})-cv(thirdpt(i,j)));
   end if;
  end if;
 end function;
 sc:=[ [ Eltseq(scv(i,j)) : j in D ] : i in D ];
// inner product in the Matsuo algebra
 ip:=function(x,y)
  a:=Eltseq(x); 
  b:=Eltseq(y);
  return 1/2*c*(&+[ a[i]*b[i] : i in D ])+
   1/8*c*h*(&+[ a[i]*b[j]*MMt[i,j] : i,j in D | j ne i ]);
 end function;
 return  Algebra< R, #D | sc >,ip;
end function;