Application Center - Maplesoft

App Preview:

Difforms2

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application




# Last edit: Nov 2 2010

difforms2 := module()
  option package;  
  description `Extension of the package difforms`;

  export
    About,
    iota, sortform, formbasis, coeffeqn,
    image, pullback, infpullback, imaginaryfromreal,
    flat, sharp, omegaflat, omegasharp,
    metriconforms, hodge, hodgeodd, omegaonforms, symphodge,
    AHitchin, KHitchin, LambdaHitchin,
    TwoFormToGram, GramToTwoForm,
    formIm, formRe, formconjugate,
    matrixmult,
    printsolution,printsolution2;

  local
    basischeck, formcheck, nosumds, sumd, maxdegree, mindegree;

  #####################################################################################

  About := proc()
    print(`------------------------------------------------------------------------------------`):
    print(`An extension of the Maple package "difforms"`):
    print(`Copyright 2010, by Fabian Schulte-Hengesbach`):
    print(`------------------------------------------------------------------------------------`):
  end proc:

  #######################################################################################
  formcheck := proc(alpha)
    description `Checks whether the argument is of type form, scalar or const and whether it contains any undeclared terms.`;
    local i,a;
    if not type(alpha,form) and not type(alpha,scalar) and not type(alpha,const) then
      error "The argument",alpha,"must be a differential form declared with the package difforms";
    elif not type(wdegree(alpha),integer) and not wdegree(alpha)=nonhmg then
      error "The argument",alpha,"contains undeclared terms.";
    elif wdegree(alpha)=nonhmg then
      a:=simpform(alpha);
      for i from 1 to nops(a) do
        if not type(wdegree(op(a)[i]),integer)  then
          error "The argument",op(a)[i],"contains undeclared terms.";
        end if:
      end do:
    end if:
  end proc:
  #######################################################################################
  basischeck := proc(a)
    description `Checks whether the argument is a list of one-forms.`;
    local i;
    if not type(a,list) then
      error "The argument",a,"must be a list of one-forms."
    else for i from 1 to nops(a) do
         if not type(op(a)[i],form) then
           error "The argument",a,"must be a list of one-forms."
         elif not wdegree(op(a)[i])=1 then
           error "The argument",a,"must be a list of one-forms."
         end if:
      end do:
    end if:
  end proc:
  ######################################################################################
  nosumds := proc(a)  
    description `Returns the number of summands of a p-form (after applying simpform). Solves a number of problems with nops.`;
    ######## To remember: formpart applies simpform, but scalarpart does not ##########
    formcheck(a);
    if type(a,scalar) or type(a,const) then
      return 1   
    elif wdegree(a)=nonhmg then
      return nops(simpform(a))  
    elif wdegree(formpart(simpform(a)) - op(formpart(simpform(a)))[1])=nonhmg then
      return 1   
    else
      return nops(formpart(simpform(a)));
    end if:
  end proc:
  ######################################################################################
  sumd := proc(a)
    description `Returns the summands of a form as a list. Fixes problem of op with only one summand.`;
    ## error handling in nosumds ##
    if nosumds(a) > 1 then
      return [op(a)];
    elif nosumds(a) = 1 then
      return [a];
    elif nosumds(a) = 0 then
      return [a];
    end if:
  end proc:
  ######################################################################################
  maxdegree := proc(a)
    description `Maximal degree of a summand of a non-homogeneous form.`;
    local max,i;
    formcheck(a);
    if wdegree(a)=nonhmg then
      max:=0:
      for i to nosumds(a) do
        if max < wdegree(sumd(a)[i]) then max := wdegree(sumd(a)[i]); end if:
      end do:
      return max;
    else return wdegree(a);
    end if:
  end proc:
  ######################################################################################
  mindegree := proc(a)
    description `Minimal degree of a summand of a non-homogeneous form.`;
    local min,i;
    formcheck(a);
    if wdegree(a)=nonhmg then
      min:=100:
      for i to nosumds(a) do
        if min > wdegree(sumd(a)[i]) then min := wdegree(sumd(a)[i]); end if:
      end do:
      return min;
    else return wdegree(a);
    end if:
  end proc:
  ######################################################################################
  iota := proc(alpha,beta)
    description `Contraction of a p-form beta with a k-vector alpha, k <= p. Both arguments have to be expressed with respect to a basis of one-forms.`;
    options remember;
    local match, output, i, sign, a, b;

    formcheck(alpha);formcheck(beta);
    a:=simpform(alpha);b:=simpform(beta);
    if alpha = 0 or beta = 0 then
      return 0
    elif maxdegree(a) > mindegree(b) then
       error "The degree of the first argument",alpha,"has to be less or equal than that of the second argument",beta;
    ### Case 0: The first argument is a scalar or a constant. ###
    elif type(a,const) or type(a,scalar) then
      return a * b;
    ### Case 1: The first argument has multiple summands -> recursion ###
    elif nosumds(a) > 1 then
       output := 0;
       for i from 1 to nosumds(a) do
          output := output + iota(sumd(a)[i],b)
       end do:
       return simpform(output);
    ### Case 2: The second argument has multiple summands -> recursion ###
    elif nosumds(b) > 1 then
         output := 0;
         for i to nosumds(b) do
            output := output+iota(a,sumd(b)[i])
         end do:
       return simpform(output);
    ### Case 3: The first argument has degree greater than 1 -> recursion ###
    elif wdegree(a) > 1 then
       if nops(formpart(a)) < 2 then
           error "The first argument",alpha,"must be expressed in a basis of one-forms."
       else
          output := scalarpart(a) * b;
          for i from 1 to nops(formpart(a)) while not output = 0 do
            output := iota(op(formpart(a))[i],output);
          end do:
         return simpform(output);
       end if:
    ### Case 4: Both arguments a and b consists of one summand and wdegree(a)=1 ###
    else
        match:=false;
        sign:=1;
        output := 1;
        if nops(formpart(b)) = 1 and wdegree(b) > 1 then
            error "The second argument",beta,"must be expressed in a basis of one-forms.";
          elif nops(formpart(b)) = 1 and formpart(b) = formpart(a) then
            return scalarpart(b)*scalarpart(a);
          elif nops(formpart(b)) = 1 then
             return 0;
          else
             for i to nops(formpart(b)) do
                if op(formpart(b))[i] = formpart(a) then
                   match := true;
                   sign := (-1)^(i+1);
                else
                   output := &^(output, op(formpart(b))[i]);
                end if:
             end do;
             if match then
                  return sign*scalarpart(b)*scalarpart(a)*output;
               else
                  return 0;
             end if:
        end if:
     end if:

   end proc:
   ######################################################################################
   sortform := proc(alpha,ob)
     description `The argument alpha must be a p-form expressed w.r.t. the ordered basis of one-forms passed in the second argument. Sortform rearranges the one-forms of each summand of alpha in the order of the ordered basis. Additionally, sortform applies simplify to the scalarpart of each summand.`;

     local i,a,b,newform,swaps,swapped,gap,tmp;

     formcheck(alpha); basischeck(ob);
     a:=simpform(simplify(alpha));
     ### Case 0: Degree 0 ##     
     if type(a,scalar) or type(a,const) then
            return a;   
     ### Case 1: Multiple summands -> recursion ##
     elif nosumds(a) > 1 then
       defform(newform=-1);
       newform:=0;
       for i from 1 to nosumds(a) do
         newform := newform + sortform(op(a)[i],ob);
       end do:
       return newform;
     ### Case 2: One summand and degree 1 ###
     elif nosumds(a) = 1 and wdegree(a)=1 then
       if member(formpart(a),ob) then
           return a;
         else
           error "The argument",a,"must be expressed w.r.t. the basis passed in the argument",ob;
       end if:
     ### Case 3: One summand and degree > 1 ###
     else
     b:=[seq(0,i = 1..wdegree(a))];
     for i from 1 to wdegree(a) do
       if not member(op(formpart(a))[i],ob,'pos') then
           error "The argument",a,"must be expressed w.r.t. the basis passed in the argument",ob;
       else b[i]:=pos;
       end if:
     end do:
     
     ### CombSort of the list b (need to count swaps for the sign) is improved bubble sort ###
     ### seems to be a fast algorithm which sorts by exchanging and uses no extra memory ###
     swaps := 0;
     gap := nops(b);
     swapped := false;
     while not gap = 1 or swapped do  
         if gap > 1 then
             gap := trunc(gap / 1.3);
             if gap = 10 or gap = 9 then
                 gap := 11;
             end if:
         end if:
         swapped := false;
         for i from 1 while i-1 + gap < nops(b) do
             if b[i] > b[i+gap] then
                 tmp:=b[i]; b[i]:=b[i+gap]; b[i+gap]:=tmp;  ##swap(b[i], b[i+gap])
                 swapped := true;
              swaps:=swaps + gap + gap-1;   ## Number of swaps
             end if:
         end do:
     end do:

     newform:=1;
     for i from 1 to wdegree(a) do
        newform:=&^(newform,ob[b[i]]);
     end do;
       return simplify(scalarpart(a)) * (-1)^swaps * newform;
     end if:
  end proc:
  ######################################################################################
  AHitchin := proc(a,ob)
    description `The procedure AHitchin assigns to a p-form alpha the (n-p)-form A(alpha) such that A(alpha) hook ob[1] wedge ... wedge ob[n] = alpha. In fact, it is A(alpha) = (-1)^p(n-p) alpha hook (ob[1])* wedge ... wedge (ob[n])*.`;
    local i, vol;

    formcheck(a); basischeck(ob);
    if wdegree(a) = nonhmg then
      error "The argument",a,"must be homogeneous.";
    elif wdegree(a) > nops(ob) then
      error "The argument",a,"must be expressed w.r.t. the basis passed in the argument",ob;
    else
    vol:=1;
    for i from 1 to nops(ob) do
      vol:=&^(vol,ob[i]);
    end do:
    return (-1)^(wdegree(a)*(nops(ob)-wdegree(a)))*iota(a,vol);
    end if:
  end proc:
  ######################################################################################
  KHitchin := proc(psi,ob)
    description `KHitchin Assigns to a m-form psi in dimension 2m the matrix representing the endomorphism X mapsto A(X hook psi wedge psi) w.r.t. to the basis passed in the second argument. Lambda^n V^* is identified with R with the help of the volume form ob[1] wedge ... wedge ob[2m].`;
    local i,j,n,output,tmp;

    formcheck(psi); basischeck(ob);
    n:=nops(ob);
    if not 2*wdegree(psi) = n then
      error "The argument", psi, "has to be m-form and the length of the basis",ob,"must be 2m.";
    else
     output:=Matrix(n);  
     for i from 1 to n do
       tmp:=AHitchin(&^(iota(ob[i],psi),psi),ob);
       if not tmp = 0 then
         for j from 1 to nosumds(tmp) do
             if member(formpart(sumd(tmp)[j]),ob,'pos') then
               output[pos,i]:=scalarpart(simpform(sumd(tmp)[j]));
             else
               error "Should never have happened"
             end if:
         end do:           
        end if:
      end do:
    end if:
    return output;
 
  end proc:
  ######################################################################################
  LambdaHitchin:=proc(psi,ob)
     description `Returns 1/m*trace(KHitchin(psi)^2) for a m-form psi in dimension 2m as a real number where Lambda^2m is identified with R with the help of the volume form ob[1] wedge ... wedge ob[2m].`;
     local i, j, K, tr;
     ### Type check is left to KHitchin ###
     ### The trace of K^2 is explicitly implemented due to inconsistencies between LinearAlgebra and difforms ###
     K:=KHitchin(psi,ob);
     tr:=0;
     for i from 1 to nops(ob) do
       for j from 1 to nops(ob) do
          tr:=tr + K[i,j]*K[j,i];
        end do:
     end do:
     
     return (1/nops(ob))*factor(tr);
  end proc:
  #######################################################################################
  TwoFormToGram:=proc(omega,ob)
     local output,n,i,j,tmp;
     description `Returns the Gram matrix of a given two-form with respect to an ordered base.`;
     
     formcheck(omega); basischeck(ob);
     if not wdegree(omega) = 2 then
        error "The first argument has to be a two-form."
     else
       n:=nops(ob);
       tmp:=sortform(omega,ob); # Checks whether omega is expressed w.r.t. ob
       output:=Matrix(n);
       for i from 1 to nosumds(omega) do
          member(op(formpart( sumd(tmp)[i] ))[1],ob,'pos1');
          member(op(formpart(sumd(tmp)[i]))[2],ob,'pos2');
          output[pos1,pos2]:=scalarpart(simpform(sumd(tmp)[i]));
          output[pos2,pos1]:=-scalarpart(simpform(sumd(tmp)[i]));
       end do:
     return output;
     end if:

  end proc:
 ######################################################################################
 GramToTwoForm:=proc(Omega,ob)
     local output,n,i,j,tmp;
     description `Converts a skew-symmetric matrix into a two-form with respect to the given ordered basis of one-forms.`;
     
     ### type check ###
     if not type(Omega,Matrix) then
        error "The first argument has to be a skewsymmetric matrix."
     elif not IsMatrixShape(Omega, skewsymmetric) then
        error "The first argument has to be a skewsymmetric matrix."
     elif not type(ob,list) then
        error "The second argument has to be a list of one-forms."
     else
        n:=nops(ob);  ## n = dimension
        for i from 1 to n do
           if not type(op(ob)[i],form) then
             error "The second argument has to be a list of one-forms."
           elif not wdegree(op(ob)[i])=1 then
             error "The second argument has to be a list of one-forms."
           end if:
        end do:
        if not LinearAlgebra[Dimension](Omega)=(n,n) then
           error "The matrix has to be of the same size as the basis."
        else
           for i from 1 to n do for j from 1 to n do
              if not type(Omega[i,j],const) and not type(Omega[i,j],scalar) then
                error "The entries of the matrix have to be zero-forms or constants,
                  i.e. type scalar or type const in the package difforms";
              end if;
           end do: end do:
        end if:
     end if:
     ### end type check ###

     defform(output=-1); output := 0;
     for i from 1 to n do
       for j from 1 to n do
          output := output + simplify((1/2)*Omega[i,j]) * &^(ob[i],ob[j]);
       end do;
     end do;
    
     return sortform(output,ob);

  end proc:
 ######################################################################################
 image := proc(v,M,ob)
    description `The Matrix M is applied to a vector v which is given as a linear combination in the base "ob" and also returned as a linear combination in ob.`;
    local i,j,n,output,pos;
    option remember;

    basischeck(ob);n:=nops(ob);
    if not type(v,form) then
        error "The first argument has to be a one-form."
    elif not wdegree(v)=1 then
        error "The first argument has to be a one-form."
    elif not type(M,Matrix) then
        error "The second argument has to be a matrix."
    elif not LinearAlgebra[Dimension](M)=(n,n) then
        error "The matrix has to be of the same size as the basis."
    end if:
    output:=0;
    if nosumds(v) > 1 then
      for i from 1 to nosumds(v) do
        output := output + image( sumd(v)[i] , M, ob);  # Recursion
      end do:
      return sortform(output,ob):
    else
      if member(formpart(v),ob, 'pos') then
         for j from 1 to n do
            output:=output + simplify(M[j,pos]) * ob[j];
         end do:
      else
        error "The first argument has to be a one-form expressed in the basis given in the third argument.";
      end if:
      return scalarpart(simpform(v)) * output;
    end if:

  end proc:
######################################################################################
pullback := proc(alpha,M,ob)
   description `Given a p-form alpha, a matrix M and a basis ob such that M represents an endomorphism w.r.t. the basis ob, the pullback of the p-form by the endomorphim is returned.`;
   local i,n,output,tmp,a,S;
   option remember;
   
   formcheck(alpha); basischeck(ob); n:=nops(ob);
   if not type(M,Matrix) then
     error "The second argument has to be a Matrix."
   elif not LinearAlgebra[Dimension](M)=(n,n) then
     error "The Matrix has to be of the same size as the basis."
   elif wdegree(alpha) = 0 then
     return alpha;
   end if:
   output:=0; tmp := 1; a:=sortform(alpha,ob);
   S:=LinearAlgebra[Transpose](M); # Act with transpose on one-forms
   ### Case 1: Multiple summands ###
   if nosumds(a) > 1 then
      for i from 1 to nosumds(a) do
        output := output + pullback(sumd(a)[i], M, ob);  # reduce to single summand
      end do:
      return sortform(output,ob):
  ### Case 2: One summand and degree = 1 ###
  elif nosumds(a) = 1 and wdegree(a)=1 then
     if member(formpart(a),ob) then
       return image(a,S,ob);
     else
       error "The argument",alpha,"must be expressed w.r.t. the basis passed in the argument",ob;
     end if:
  ### Case 3: One summand and degree > 1 ###
  else
     for i from 1 to wdegree(a) do
       if not member(op(formpart(a))[i],ob,'pos') then
           error "The argument",alpha,"must be expressed w.r.t. the basis passed in the argument",ob;
       else tmp:=sortform(&^(tmp,image(ob[pos],S,ob)),ob);
       end if:
     end do:
     return scalarpart(a)*tmp;
   end if:

end proc:
#########################################################################
infpullback := proc(alpha,M,ob)
   description `Given a p-form alpha, a basis ob and a matrix M representing an endomorphism w.r.t. this basis, the action of the endomorphism as a derivation on the p-form is computed.`;
   local i,n,output,tmp1,tmp2,a,S;
   option remember;
   
   formcheck(alpha); basischeck(ob); n:=nops(ob);
   if not type(M,Matrix) then
     error "The second argument has to be a Matrix."
   elif not LinearAlgebra[Dimension](M)=(n,n) then
     error "The Matrix has to be of the same size as the basis."
   elif wdegree(alpha) = 0 then
     return alpha;
   end if:
   output:=0; a:=sortform(alpha,ob);
   S:=LinearAlgebra[Transpose](M); # Act with transpose on one-forms
   ### Case 1: Multiple summands ###
   if nosumds(a) > 1 then
      for i from 1 to nosumds(a) do
        output := output + infpullback(sumd(a)[i], M, ob);  # reduce to single summand
      end do:
      return sortform(output,ob):
  ### Case 2: One summand and degree = 1 ###
  elif nosumds(a) = 1 and wdegree(a)=1 then
     if member(formpart(a),ob) then
       return image(a,S,ob);
     else
       error "The argument",alpha,"must be expressed w.r.t. the basis passed in the argument",ob;
     end if:
  ### Case 3: One summand and degree > 1 ###
  else
     for i from 1 to wdegree(a) do
       if not member(op(formpart(a))[i],ob,'pos') then
           error "The argument",alpha,"must be expressed w.r.t. the basis passed in the argument",ob;
       else
         tmp1:=iota(ob[pos],formpart(a));
         tmp2:=sortform(&^(image(ob[pos],S,ob),tmp1),ob);
         output:= output + scalarpart(a)*tmp2;
       end if:
     end do:
     return sortform(output,ob);
   end if:

end proc:
############################################################
imaginaryfromreal := proc(argument_a,J,orderedbase)
   local i,j,ausgabe,n,k,l,coeff,baseform,psi;
   description `Computes the imaginary part of a given three-form in dimension 6 with lambdaHitchin < 0.`;
   ### Evaluates the identity iota_X(imaginarypart)=-iota_JX(realpart) on the induced basis of three-forms ###
   ### J could be computed from argument_a in the procedure. Since it may look VERY complicated the user is forced  to compute in before with the proceudre KHitchin.  ###
   if not wdegree(argument_a) = 3 then error "The first argument has to be a three-form."; end if:
   defform(baseform=-1,ausgabe=-1,psi=-1); ausgabe:=0; n:=6; baseform := 1;
   psi:=sortform(argument_a,orderedbase);
   for i from 1 to n do
     for j from i+1 to n do
       for k from j+1 to n do             
          ### Koeffizient von e_i,e_j,e_k im Imaginaerteil: (Je_i,e_j,e_k) einsetzen in Realteil und Minus ###
          coeff:=simplify(iota(&^(image(orderedbase[i], J, orderedbase), orderedbase[j], orderedbase[k]), psi));
          if not coeff=0 then
            baseform := &^ (orderedbase[i],orderedbase[j],orderedbase[k]);
            ausgabe := ausgabe - coeff * baseform; ## Der Koeffizient MIT Minuszeichen muß also vor e_i,e_j,e_k stehen ##
          end if:
       end do:
     end do:
   end do:
   return sortform(ausgabe,orderedbase);
end proc:
######################################################################################
formRe := proc(alpha)
  local newform,a,j;
       description `Computes the real part of a complex k-form.`;
  ## type check in nosumds ##
  defform (newform=0,a=0); a := simpform(alpha); newform := 0;
  
  for j from 1 to nosumds(a) do
     if not Re(scalarpart(simpform(sumd(a)[j]))) = 0 then
       newform := newform + Re(scalarpart(simpform(sumd(a)[j])))*formpart(sumd(a)[j]);
     end if:
  end do:
  return simpform(newform);
end proc:
######################################################################################
formIm := proc(alpha)
  local newform,a,j;
       description `Computes the imaginary part of a complex k-form.`;
  ## type check in nosumds ##
  defform (newform=0,a=0); a := simpform(alpha); newform := 0;
  
  for j from 1 to nosumds(a) do
     if not Im(scalarpart(simpform(sumd(a)[j]))) = 0 then
       newform := newform + Im(scalarpart(simpform(sumd(a)[j])))*formpart(sumd(a)[j]);
     end if:
  end do:
return simpform(newform);
end proc:
######################################################################################
formconjugate := proc(alpha)
  description `Computes the complex conjugate of a complex k-form.`;
  ## type check in formRe and formIm ##
  return formRe(alpha) - simpform(I*formIm(alpha));
end proc:
######################################################################################
coeffeqn := proc(form_a,form_b,ob)
     local i,j,a,b,n_a,n_b,stilltocheck,found,counter,diff,arethesame,returnset;
     description `Compares the coefficients of two p-Forms: Returns the set of the simplified coefficients of the difference of the two forms.`;

     ## type check in sortform ##
     defform(a=0,b=0);
     a:=sortform(form_a,ob);b:=sortform(form_b,ob);
     n_a:=nosumds(a);n_b:=nosumds(b);
     stilltocheck:={seq(i,i=1..n_b)};
     counter:=1;
     arethesame:=true;
     returnset:={};
   
     for i from 1 to n_a do
        found:=false;
        for j in stilltocheck do
          if formpart(sumd(a)[i])=formpart(sumd(b)[j]) then
            diff:=simplify(scalarpart(simpform(sumd(a)[i]))-scalarpart(simpform(sumd(b)[j])));
            if not diff = 0 then
              arethesame:=false;
              returnset:=returnset union {diff=0};
            end if:
            found:=true;
            stilltocheck:=stilltocheck minus {j};
            counter:=counter+1:
          end if:
        end do:
        if not found then
            returnset:=returnset union {simplify(scalarpart(simpform(sumd(a)[i])))=0};
            counter:=counter+1:
            arethesame:=false;
        end if:
     end do:
     for i in stilltocheck do
            returnset:=returnset union {simplify(scalarpart(simpform(sumd(b)[i])))=0};
            counter:=counter+1:
            arethesame:=false;
     end do:
     
     if arethesame then
         print(`The forms are identical!`);
       else
         #print(nops(returnset),`Non-trivial equations:`);
         return returnset;
     end if:
  end proc:

  #####################################################
  omegaonforms := proc(alpha,beta,omega,ob)
    ### Remark: Performance could easily be improved, e.g. the Matrix Omega is inverted with every recursion again ###
    description `Computes the pairing of the p-forms alpha and beta induced by an (almost) symplectic form omega. All three forms have to be expressed w.r.t. to the basis of one-forms ob.`;
    local a,b,i,j,tmp,Omega,Odual,odual, output;
    option remember;
    
    a:=sortform(alpha,ob); b:=sortform(beta,ob); sortform(omega,ob);
    ## sortform contains all necessary type checks except: ##
    if wdegree(alpha) = nonhmg or wdegree(beta) = nonhmg then
       error "The arguments", alpha, "and", beta,"must be homogeneous."
    elif alpha = 0 or beta = 0 then
       return 0;
    elif not wdegree(alpha) = wdegree(beta) then
       error "The arguments", alpha, "and", beta,"must be of the same degree."
    elif not wdegree(omega) = 2 then
       error "The argument",omega,"must be a two-form."
    end if:
    Omega:=TwoFormToGram(omega,ob);
    if LinearAlgebra[Determinant](Omega)=0 then
       error "The argument",omega,"must be non-degenerate."
    ### Case 0: Both arguments have degree 0 ####
    elif wdegree(alpha) = 0 and wdegree(beta) = 0 then
       return simplify(alpha * beta);
    end if:
    Odual:=LinearAlgebra[Transpose](LinearAlgebra[MatrixInverse](Omega));
    odual := GramToTwoForm(Odual, ob); # dual 2-vector of omega, i.e. induced symp. form on one-forms
    ### Case 1: Both arguments have degree 1 ###
    if wdegree(alpha) = 1 and wdegree(beta) = 1 then
      return(iota(&^(alpha,beta),odual));
    ### Case 2: Argument alpha has multiple summands ###
    elif nosumds(a) > 1 then
      output:=0:
      for i from 1 to nosumds(a) do
        output := output + omegaonforms(sumd(a)[i], b, omega, ob);  # reduce to single summand
      end do:
      return sortform(output,ob):
    ### Case 3: Argument beta has multiple summands ###
    elif nosumds(b) > 1 then
      output:=0:
      for i from 1 to nosumds(b) do
        output := output + omegaonforms(a,sumd(b)[i], omega, ob);  # reduce to single summand
      end do:
      return sortform(output,ob):
    ### Case 4: Both arguments are decomposable, i.e. consist of a single summand expressed w.r.t. ob ###
    else # apply definition of the induced pairing on decomposable forms
      tmp:=Matrix(wdegree(a));
      for i from 1 to wdegree(a) do for j from 1 to wdegree(a) do
        tmp[i,j]:=iota(&^(op(formpart(a))[i],op(formpart(b))[j]),odual);
      end do: end do:
      return simplify(scalarpart(a)*scalarpart(b)*LinearAlgebra[Determinant](tmp));
    end if:
end proc:
#########################################################################################
formbasis := proc(k,ob)
  description `Returns the basis of k-forms defined by the basis of one-forms ob as a list (sorted lexicographically).`;
  local i,j,n,basis,index,tmp:
  basischeck(ob);n:=nops(ob);
  if not type(k,nonnegint) then
    error "The argument",k,"must be a non-negative integer."
  end if:
  index:=combinat[choose](n,k); # subsets of { 1.. n } of order k
  basis:=[]:
  for i from 1 to nops(index) do
    tmp:=1:
    for j from 1 to k do
      tmp:=&^(tmp,ob[index[i,j]]);
    end do:
    basis:=[op(basis),tmp]
  end do:
  return basis;
end proc:
#########################################################################################
symphodge := proc(alpha,omega,ob)
  description `Returns the image of the k-form alpha under the symplectic hodge operator defined by the symplectic form omega. Both alpha and omega must be expressed w.r.t. the basis of 1-forms ob.`;
  local a,k,i,n,basis,vol,output:
  option remember;

  a:=sortform(alpha,ob); sortform(omega,ob); # includes type-check
  basischeck(ob); n:=nops(ob);
  if not wdegree(omega) = 2 then
    error "The argument",omega,"must be a two-form."
  elif not n mod 2 = 0 then
    error "The dimension must be even."
  elif wdegree(a)=nonhmg then
    error "The argument",alpha,"must be homogeneous.";
  elif not type(wdegree(a),nonnegint) then
    error "Should not happen.";
  end if:
  k:=wdegree(alpha); vol:=1/(factorial(n/2));
  for i from 1 to n/2 do
    vol:=vol &^ omega;
  end do:
  vol:=sortform(vol,ob);
  basis:=formbasis(k,ob);
  output:=0:
  for i from 1 to nops(basis) do
    output := output + omegaonforms(basis[i],a,omega,ob)*iota(basis[i],vol)
  end do:
  return sortform(output,ob);
end proc:
#########################################################################################
 metriconforms := proc(alpha,beta,Gdual,ob)
    description `Computes the inner product of the p-forms alpha and beta induced by the inner product g which is passed as a symmetric Matrix Gdual representing g w.r.t. to the basis of one-forms ob. The arguments alpha and beta be expressed w.r.t. to the basis of one-forms ob.`;
    local a,b,i,j,n,tmp,G,output,pos1,pos2;
    option remember;
    
    a:=sortform(alpha,ob); b:=sortform(beta,ob);
    ## sortform contains all necessary type checks except: ##
    G:=simplify(Gdual); n:=nops(ob);
    if wdegree(alpha) = nonhmg or wdegree(beta) = nonhmg then
       error "The arguments", alpha, "and", beta,"must be homogeneous."
    elif not type(G,Matrix) or not IsMatrixShape(G, symmetric) then
       error "The third argument must be a symmetric Matrix."
    elif not LinearAlgebra[Dimension](G)=(n,n) then
       error "The Gram Matrix of the metric has to be of the same size as the basis."
    elif LinearAlgebra[Determinant](G)=0 then
       error "The metric must be non-degenerate."
    elif alpha = 0 or beta = 0 then
       return 0;
    elif not wdegree(alpha) = wdegree(beta) then
       error "The arguments", alpha, "and", beta,"must be of the same degree."
    end if:

    ### Case 0: Both arguments have degree 0 ####
    if wdegree(alpha) = 0 and wdegree(beta) = 0 then
       return simplify(alpha * beta);
    end if:
    ### Case 1: Argument alpha has multiple summands ###
    if nosumds(a) > 1 then
      output:=0:
      for i from 1 to nosumds(a) do
        output := output + metriconforms(sumd(a)[i], b, G, ob);  # reduce to single summand
      end do:
      return sortform(output,ob):
    ### Case 2: Argument beta has multiple summands ###
    elif nosumds(b) > 1 then
      output:=0:
      for i from 1 to nosumds(b) do
        output := output + metriconforms(a,sumd(b)[i], G, ob);  # reduce to single summand
      end do:
      return sortform(output,ob):
    ### Case 1: Both arguments have degree 1 (and consist of a single summand.) ###
    elif wdegree(a) = 1 and wdegree(b) = 1 then
      member(formpart(a),ob,'pos1');
      member(formpart(b),ob,'pos2');
      return(simplify(scalarpart(a)*scalarpart(b)*G[pos1,pos2]));
    ### Case 3: Both arguments are decomposable and of degree > 1, ###
    ### i.e. consist of a single summand expressed w.r.t. ob ###
    else # apply definition of the induced pairing on decomposable forms
      tmp:=Matrix(wdegree(a));
      for i from 1 to wdegree(a) do
        member(op(formpart(a))[i],ob,'pos1');
        for j from 1 to wdegree(a) do
          member(op(formpart(b))[j],ob,'pos2');
          tmp[i,j]:=G[pos1,pos2];
      end do: end do:
      return simplify(scalarpart(a)*scalarpart(b)*LinearAlgebra[Determinant](tmp));
    end if:
end proc:
#################################################################################
hodge := proc(alpha,Gdual,ob)
  description `Returns the image of the k-form alpha under the hodge operator defined by the inner product g which is passed as a symmetric Matrix Gdual representing g w.r.t. to the oriented basis of one-forms ob. The argument alpha must also be expressed w.r.t. the basis of 1-forms ob. WARNING: Due to incompatibilities between absolute value and the package difforms, the implemented procedure only returns the correct signs if the number of minuses in the signature of g is even, i.e. if det(G)>0. Otherwise, use "hodgeodd".`;
  local a,k,i,n,G,basis,vol,output:
  option remember;

  a:=sortform(alpha,ob); # includes type-checks
  n:=nops(ob);
  G:=simplify(Gdual);
  if not type(G,Matrix) or not IsMatrixShape(G, symmetric) then
       error "The second argument must be a symmetric Matrix."
  elif LinearAlgebra[Determinant](G)=0 then
       error "The metric must be non-degenerate."
  elif wdegree(a)=nonhmg then
    error "The argument",alpha,"must be homogeneous.";
  elif not type(wdegree(a),nonnegint) then
    error "Should not have happened after type-checks.";
  end if:
  k:=wdegree(alpha);
  vol := simplify(sqrt(1/LinearAlgebra[Determinant](G)))*formbasis(n, ob)[1];
  basis:=formbasis(k,ob);
  output:=0:
  for i from 1 to nops(basis) do
    output := output + metriconforms(basis[i],a,G,ob)*iota(basis[i],vol)
  end do:
  return sortform(output,ob);
end proc:
#########################################################################################
#################################################################################
hodgeodd := proc(alpha,Gdual,ob)
  description `Returns the correct hodge dual if the number of minuses in the signature of g is odd, i.e. if det(G)<0. `;
  local a,k,i,n,G,basis,vol,output:
  option remember;

  a:=sortform(alpha,ob); # includes type-checks
  n:=nops(ob);
  G:=simplify(Gdual);
  if not type(G,Matrix) or not IsMatrixShape(G, symmetric) then
       error "The second argument must be a symmetric Matrix."
  elif LinearAlgebra[Determinant](G)=0 then
       error "The metric must be non-degenerate."
  elif wdegree(a)=nonhmg then
    error "The argument",alpha,"must be homogeneous.";
  elif not type(wdegree(a),nonnegint) then
    error "Should not have happened after type-checks.";
  end if:
  k:=wdegree(alpha);
  vol := simplify(sqrt(-1/LinearAlgebra[Determinant](G)))*formbasis(n, ob)[1];
  basis:=formbasis(k,ob);
  output:=0:
  for i from 1 to nops(basis) do
    output := output + metriconforms(basis[i],a,G,ob)*iota(basis[i],vol)
  end do:
  return sortform(output,ob);
end proc:
#########################################################################################
matrixmult := proc(A,B)
  description `Matrix multiplication. Necessary due to inconsistencies of the package difforms and LinearAlgebra. Arguments must be of type Matrix or const or scalar.`;
  local i,j,k,n,output:
  if not type(A,Matrix) and not type(A,scalar) and not type(A,const) then
      error "The argument",A,"must be of type Matrix or const or scalar";
  elif not type(B,Matrix) and not type(B,scalar) and not type(B,const) then
      error "The argument",B,"must be of type Matrix or const or scalar";
  elif type(B,Matrix) and type(A,scalar) or type(A,const) then
    output:=Matrix(LinearAlgebra[RowDimension](B),LinearAlgebra[ColumnDimension](B));
    for i from 1 to LinearAlgebra[RowDimension](B) do
      for j from 1 to LinearAlgebra[ColumnDimension](B) do
        output[i,j]:=simplify(A*B[i,j]);
      end do:
    end do:
    return(output);
  ### Rightmultiplication with scalar: Assumes commutatitivy of ground field ###
  elif type(A,Matrix) and type(B,scalar) or type(B,const) then
    return(matrixmult(B,A));
  elif type(A,Matrix) and type(B,Matrix) then
     if not LinearAlgebra[ColumnDimension](A) = LinearAlgebra[RowDimension](B) then
        error "The number of colums of the first argument must match the number of rows of the second argument.";
     end if:
     output:=Matrix(LinearAlgebra[RowDimension](A),LinearAlgebra[ColumnDimension](B)):
     n:=LinearAlgebra[RowDimension](B):
     for i from 1 to LinearAlgebra[RowDimension](A) do
      for j from 1 to LinearAlgebra[ColumnDimension](B) do
        output[i,j]:=0:
        for k from 1 to n do
          output[i,j]:=output[i,j]+A[i,k]*B[k,j]:
        end do:
      end do:
    end do:
     return(output);
    end if:
end proc:
##########################################################################################
flat := proc(v,G,ob)
  description `Returns the one-form which is the metric dual of the vector v w.r.t to the metric G.`;
  local S;
  formcheck(v);
  if wdegree(v)=1 then
    return image(v,G,ob);
  else
    S:=LinearAlgebra[Transpose](G); #Use "pullback" with transpose to push forward
    return pullback(v,S,ob);
  end if:
end proc:
##########################################################################################
sharp := proc(alpha,Gdual,ob)
  description `Returns the vector which is the metric dual of the one-form alpha w.r.t to the metric g which is passed in the basis of one-forms ob.`;
  local S;
  formcheck(alpha);
  if wdegree(alpha)=1 then
    return image(alpha,Gdual,ob);
  else
    S:=LinearAlgebra[Transpose](Gdual); #Use "pullback" with transpose to push forward
    return pullback(alpha,S,ob);
  end if:
end proc:
##########################################################################################
omegaflat := proc(v,omega,ob)
  description `Returns the one-form which is the symplectic dual of the vector v w.r.t to the symplectic form omega.`;
  local O,S;
  O:=LinearAlgebra[Transpose](TwoFormToGram(omega,ob)); # Transpose due to convection \flat(X)(Y)=\omega(X,Y)
  formcheck(v);
  if wdegree(v)=1 then
    return image(v,O,ob);
  else
    S:=LinearAlgebra[Transpose](O); #Use "pullback" with transpose to push forward
    return pullback(v,S,ob);
  end if:
end proc:
##########################################################################################
omegasharp := proc(alpha,omega,ob)
  description `Returns the vector which is the metric dual of the one-form alpha w.r.t to the metric g which is passed in the basis of one-forms ob.`;
  local Oinverse,S;
  formcheck(alpha);
  Oinverse:=LinearAlgebra[Transpose](LinearAlgebra[MatrixInverse](TwoFormToGram(omega,ob)));
  if wdegree(alpha)=1 then
    return image(alpha,Oinverse,ob);
  else
    S:=LinearAlgebra[Transpose](Oinverse); #Use "pullback" with transpose to push forward
    return pullback(alpha,S,ob);
  end if:
end proc:
##########################################################################################
printsolution2 := proc(sol)
   local i,returnstring,anotherstring,txt,txt2,x;
   description `Prepares a solution for copy and paste including unassign.`;
   ### type check ###
   if not type(sol,set) then
     error `The argument has to be a set of equations which will be turned into assignments.`;
   end if:
   returnstring:="";anotherstring:="unassign(";
   for i from 1 to nops(sol) do
     txt:= convert(op(sol)[i],string);
     x:=SearchText("=", txt);
     if not substring(txt, 1 .. x-2) = substring(txt, x+2 .. length(txt)) then
        if anotherstring <> "unassign(" then anotherstring:=cat(anotherstring,", "); end if;
        anotherstring:=cat(anotherstring,"'",substring(txt,1..x-2),"'");
        txt:= StringTools[Delete](txt,x-1..x+1);
        txt := StringTools[Insert](txt, x-2, ":=");
        txt := StringTools[Insert](txt, length(txt), ": ");
        returnstring:=cat(returnstring,txt);
     end if:
   end do:
  printf("%s\n",returnstring);
  anotherstring:=cat(anotherstring,");");
  printf("%s\n",anotherstring);
 end proc:
################################################################################################
printsolution := proc(sol)
   local i,returnstring,txt,txt2,x;
   description `Prepares a solution for copy and paste.`;
   ### type check ###
   if not type(sol,set) then
     error `The argument has to be a set of equations which will be turned into assignments.`;
   end if:
   returnstring:="";
   for i from 1 to nops(sol) do
     txt:= convert(op(sol)[i],string);
     if not substring(txt, 1 .. SearchText("=", txt)-2) = substring(txt, SearchText("=", txt)+2 .. length(txt)) then
        x := SearchText("=", txt);
        txt := StringTools[Delete](txt,x-1..x+1);
        txt := StringTools[Insert](txt, x-2, ":=");
        txt := StringTools[Insert](txt, length(txt), ": ");
        returnstring:=cat(returnstring,txt);
     end if:
   end do:
  printf("%s",returnstring);
 end proc:

end module:

#Create library once with "march" in a folder of your choice

# Add the line -- libname:="e:\\home\\maplelib",libname: -- to maple.ini