// Magma program script RealQdrPxP.m, 2015/05/14

intrinsic RealQdrPxP(d::RngIntElt,p::RngIntElt){}

   SetClassGroupBounds("GRH");

   a := false; // admissible
   if (1 eq d mod 4) and (1 lt d) and IsSquarefree(d) then
      a := true;
   end if; // 1 eq d
   if (0 eq d mod 4) then
      r := d div 4; // residue
      if ((2 eq r mod 4) or (3 eq r mod 4)) and IsSquarefree(r) then
         a := true;
      end if; // r mod 4
   end if; // 0 eq d

   if (true eq a) then
      ZX<X> := PolynomialRing(Integers());
      K := NumberField(X^2-d); // base field
      // K := QuadraticField(d);
      O := MaximalOrder(K);
      C,mC := ClassGroup(O);

      if ([p,p] eq pPrimaryInvariants(C,p)) then
         printf "%7o: ",d;

         sS := Subgroups(C: Quot := [p]);
         sI := [];
         for j in [1..p+1] do
            Append(~sI,0);
         end for; // j
         n := Ngens(C);
         ct := 0; // local counter
         for x in sS do
            ct := ct+1;
            if (Order(C.(n-1)) div p)*C.(n-1) in x`subgroup then
               sI[1] := ct;
            end if; // n-1
            if (Order(C.n) div p)*C.n in x`subgroup then
               sI[2] := ct;
            end if; // n
            for e in [1..p-1] do
               if ((Order(C.(n-1)) div p)*C.(n-1))+(e*(Order(C.n) div p)*C.n) in x`subgroup then
                  sI[e+2] := ct;
               end if; // product
            end for; // e
         end for; // x

         sA := [AbelianExtension(Inverse(mQ)*mC)
                where Q,mQ := quo<C|x`subgroup>: x in sS];
         sN := [NumberField(x): x in sA];
         sR := [MaximalOrder(x): x in sA];
         sF := [AbsoluteField(x): x in sN];
         sM := [MaximalOrder(x): x in sF];
         sM := [OptimizedRepresentation(x): x in sF];
         sA := [NumberField(DefiningPolynomial(x)): x in sM];
         sO := [Simplify(LLL(MaximalOrder(x))): x in sA];

         TTT := [];
         for j in [1..#sO] do
            CO := ClassGroup(sO[j]);
            Append(~TTT,pPrimaryInvariants(CO,p));
         end for; // j

         TKT := [];
         for j in [1..#sR] do
            Collector := [];
            I := sR[j]!!mC( (Order(C.(n-1)) div p)*C.(n-1) );
            if IsPrincipal(I) then
               Append(~Collector,sI[1]);
            end if; // I
            I := sR[j]!!mC( (Order(C.n) div p)*C.n );
            if IsPrincipal(I) then
              Append(~Collector,sI[2]);
            end if; // I
            for e in [1..p-1] do
               I := sR[j]!!mC( ((Order(C.(n-1)) div p)*C.(n-1)+(e*(Order(C.n) div p)*C.n)) );
               if IsPrincipal(I) then
                  Append(~Collector,sI[e+2]);
               end if; // I
            end for; // e
            if (2 le #Collector) then
               Append(~TKT,0);
            else
               Append(~TKT,Collector[1]);
            end if; // #Collector
         end for; // j

         TAB := []; // Taussky conditions A and B
         for j in [1..#TKT] do
            if (j eq TKT[j]) or (0 eq TKT[j]) then
               Append(~TAB,"A");
            else
               Append(~TAB,"B");
            end if; // fixed point or total
         end for; // j

         printf "%o; %o; ( ",TKT,TAB;
         for j in [1..#TTT] do
            printf "%o ",TTT[j];
         end for; // j
         printf ")\n";

      end if; // type [p,p]
   end if; // admissible

end intrinsic; // RealQdrPxP