// Let V be a vector space of dimension d over GF(q) (where q is a prime power).
// Then there is a natural action of GL(d,q) on ExteriorSquare(V) and on
// DirectSum(V,ExteriorSquare(V)). Bettina Eick and I have computed the PORC
// formulae giving the numbers of orbits of k-dimensional subspaces of
// ExteriorSquare(V) and of DirectSum(V,ExteriorSquare(V)) for d = 1,2,3,4,5,6,7.
// In our paper these functions are denoted by w_{d,k}(q), w^{}_{d,k}(q).

// This Magma program gives access to our formulae.

// WedgeProduct(d,k,0) returns a string giving the PORC formula for the number of orbits 
// of k-dimensional subspaces of ExteriorSquare(V), with q symbolic. This number is the 
// number of d-generator class two Lie algebras of dimension d+k over GF(q). If q is an 
// odd prime it is also the number of d-generator class two qroups of order q^{d+k} with 
// exponent q. If you give q a value, then you can use eval to evaluate this formula at q.
// For example:

// a:=WedgeProduct(6,5,0);
// for q in [2,3,4,5] do
//   eval a;
// end for;

// If q is a prime power then WedgeProduct(d,k,q) gives the number of orbits (as an integer).

// WedgeProductPlus(d,k,0) returns a string giving the PORC formula for the number of orbits 
// of k-dimensional subspaces of DirectSum(V,ExteriorSquare(V)), with q symbolic. If q is an 
// odd prime, then this number is also the number of d-generator q-class two groups of order 
// q^{d+k}.

// If q is a prime power then WedgeProductPlus(d,k,q) gives the number of orbits (as an integer).

WP := [];
WP[1] := [];
WP[2] := [];
WP[3] := ["1"];
WP[4] := ["2","4","6"];
WP[5] := ["2","6","22","57"];
WP[5][5] := "3*q+2*GCD(q,2)+GCD(q,3)+2*GCD(q-1,3)+GCD(q-1,4)+63";

load lie6;
load lie7;

WedgeProduct := function( d,k,q )

    // these cases do not exist
    if d notin {2,3,4,5,6,7} then "d must be an integer in the range 2,3,4,5,6,7"; return 0; end if;

    // these cases do not exist
    if k notin {i:i in [1..(d*(d-1)) div 2]} then "k must be an integer in the range [1..d(d-1)/2]"; return 0; end if;

    // these cases do not exist
    if q ne 0 and not IsPrimePower(q) then "q must be 0 or a prime-power"; return 0; end if;

    // catch some trivial cases
    l := d*(d-1)/2;
    if k eq l then if q eq 0 then return "1"; else return 1; end if; end if;
    if k gt l/2 then k := l-k; end if;

    f := WP[d][k];
    if q eq 0 then return f; else; return eval f; end if;

end function;


PP := [];
PP[1] := [];
PP[2] := ["3"];
PP[3] := ["4", "q-GCD(q-0,2)+15", "3*q-3*GCD(q-0,2)+30" ];
PP[4] := [ "6", "4*q-2*GCD(q-0,2)+50" ];
PP[4][3] := 
  "q^5+2*q^4+7*q^3-q^2*GCD(q-0,2)+26*q^2-10*q*GCD(q-0,2)+q*GCD(q-1,3)"*
  "+98*q-37*GCD(q-0,2)-1/2*GCD(q-0,3)+4*GCD(q-1,3)+GCD(q-1,4)+615/2";
PP[4][4] := 
  "q^8+2*q^7+5*q^6+10*q^5+24*q^4-3*q^3*GCD(q-0,2)+56*q^3-15*q^2*GCD(q-0,2)"*
  "+q^2*GCD(q-1,3)+153*q^2-52*q*GCD(q-0,2)+6*q*GCD(q-1,3)+2*q*GCD(q-1,4)"*
  "-9/4*GCD(q-0,2)*GCD(q-1,3)-1/4*GCD(q-0,2)*GCD(q-2,3)+375*q-121*GCD(q-0,2)"*
  "-15/8*GCD(q-0,3)+171/8*GCD(q-1,3)+3/8*GCD(q-2,3)+6*GCD(q-1,4)"*
  "-1/4*GCD(q-0,5)+GCD(q-1,5)+6565/8";
PP[4][5] := 
  "q^9+2*q^8+5*q^7+10*q^6+21*q^5-q^4*GCD(q-0,2)+45*q^4-8*q^3*GCD(q-0,2)"*
  "+102*q^3-26*q^2*GCD(q-0,2)+2*q^2*GCD(q-1,3)+q^2*GCD(q-1,4)+242*q^2"*
  "-81*q*GCD(q-0,2)+10*q*GCD(q-1,3)+4*q*GCD(q-1,4)-1/4*q*GCD(q-0,5)"*
  "-4*GCD(q-0,2)*GCD(q-1,3)+2221/4*q-167*GCD(q-0,2)-GCD(q-0,3)+36*GCD(q-1,3)"*
  "+12*GCD(q-1,4)-1/4*GCD(q-0,5)+2*GCD(q-1,5)+4353/4";

load group5;
load group6;
load group7;

WedgeProductPlus := function( d,k,q )

    // these cases do not exist
    if d notin {1,2,3,4,5,6,7} then "d must be an integer in the range 1,2,3,4,5,6,7"; return 0; end if;

    // these cases do not exist
    if k notin {i:i in [1..(d*(d+1)) div 2]} then "k must be an integer in the range [1..d(d+1)/2]"; return 0; end if;

    // these cases do not exist
    if q ne 0 and not IsPrimePower(q) then "q must be 0 or a prime-power"; return 0; end if;

    // catch some trivial cases
    l := d*(d+1)/2;
    if k eq l then if q eq 0 then return "1"; else return 1; end if; end if;
    if k gt l/2 then k := l-k; end if;

    f := PP[d][k];
    if q eq 0 then return f; else return eval f; end if;

end function;


