isintegralmatrix:=function(A)
 n:=NumberOfRows(A);
 m:=NumberOfColumns(A);
 for i in [1..n] do for j in [1..m] do if not IsIntegral(A[i][j])
 then return false; end if; end for; end for;
 return true;
 end function;



 function RealSchurIndex(E)

   B:=Basis(E);
   n:=#B;
   wn := Isqrt(n) ;
   nm := wn*(wn-1) div 2;
   np := wn*(wn+1) div 2;
   F:=MatrixRing(Rationals(),n) ! 0; 
   for i in [1..n] do 
      for j in [i..n] do F[i][j]:=Trace(B[i]*B[j]);
      F[j][i]:=F[i][j];
      end for;
   end for;
   S,T:=OrthogonalizeGram(F);
   dm := 0; dp :=0; 
   for i in [1..n] do
      if S[i][i] gt 0 then dp +:=1; end if; 
      if S[i][i] lt 0 then dm +:=1; end if;
   end for;
   "n-:", nm;
   "n+:", np;
   "d-:", dm;
   "d+:", dp;
    return dp,dm;

 end function; 

function SchurIndexJac(E)

//assumes that 1 is in E //
/* 
returns the rational SchurIndex (ss) the 
sequence of local SchurIndices (one for each ramified prime) 
and the real Schur index of an algebra E
if E has uniformily distributed invariants, so the center of
E is a Galois extension K of Q and the local 
Schur index does not depend on the prime ideal of K but only
on the rational prime contained in it.
*/

     "Get Endo dim";
     time e := Dimension(E);
     "Endo dim", e;

     "Get Centre";
     time C := Centre(E);
     d := Dimension(C);
     "Centre dim", d;

     m := Isqrt(e div d);
     assert d*m^2 eq e;
     rootofdimKofE := m;

     "Get regular rep";
     time ER := RegularRepresentation(Algebra(E));
     "Get Basis";
     time B:=Basis(ER);

     is_good := func<X | Degree(MinimalPolynomial(X)) eq d>;
     alpha := 0;
     for i := 2 to d do
       if is_good(C.i) then
           alpha := C.i;
           break;
       end if;
     end for;
     if alpha eq 0 then
       range := 1;
       repeat
           alpha := &+[Random(-range, range)*X: X in Basis(C)];
           range +:= 1;
       until is_good(alpha);
     end if;

     f := MinimalPolynomial(alpha);
     K := NumberField(f);
     D := Discriminant(RingOfIntegers(K));

     /*
     R = ROI(K).
     DR = Discriminant of Mat_m(R).
     */

     DR := D^(m^2);

     "K:", K;
     "D:", D;
     "DR:", DR;

       F := MatrixRing(Integers(), e)!0;
       for i  in [1..e] do
          for j in [i..e] do
           F[i,j] := Trace(B[i]*B[j]);
           F[j,i] := F[i,j];
          end for;
       end for;
       FF := Parent(F)![ExactQuotient(x, m): x in Eltseq(F)];
     "DFF",Determinant(FF);
     DFF := ExactQuotient(Determinant(FF), DR);
     L := Factorization(DFF);

     "DFF:", DFF;
     "Factorization:", L;

realsi := 0;
//Calculate the real Schur Index, if K is totally real
// realsi = 0, if K is not totally real 
// realsi = 1, if the real Schur Index is 1
// realsi = -1, if the real Schur Index is 2
 if K eq Rationals() or IsTotallyReal(K) then
         S,T:=OrthogonalizeGram(F);
         dm:=0; dp:=0;
         for i in [1..e] do
           if S[i][i] gt 0 then dp +:=1; end if;
           // if S[i][i] lt 0 then dm +:=1; end if;
         end for;
         mm:= d*m*(m-1) div 2;
         mp:= d*m*(m+1) div 2;
          richtig := false;
         if dp eq mp then realsi := 1; richtig := true; end if;
         if dp eq mm then realsi := -1; richtig := true; end if;
         if not richtig then 
             "m :",m;
             "d :",d;
             "dp :",dp;
             "dm :",dm;
         end if;
     end if; 
     "RealSchurIndex :",realsi;
     
     if Abs(DFF) eq 1 then
       return E,realsi;
     end if;

     EO := ER;
     RM1:=B;
     RB1:=MatrixRing(Integers(),e) ! 1;


     for pt in L do
       p := pt[1];

       "Prime p:", p;
        number:=0;
        new_order := true;
        while new_order do
          number +:= 1;
        new_order:=false;
       LM1:=[]; //Left regular representation
       for i in [1..e] do X:=Parent(RM1[1]) ! 0;
                       for j,k  in [1..e] do
                       X[j][k]:=RM1[j][i][k];
               end for;
               Append(~LM1,X);
       end for;
           Ep := sub<MatrixRing(GF(p),e) | [ Matrix(GF(p),X) : X in RM1 cat LM1]>;
//take both right and left regular reps. since we are interested in twosided ideals
           Mp := RModule(Ep);
           maxsub := MaximalSubmodules(Mp);

           "#maxsub:", #maxsub;
            S:=maxsub[1];
            for si in [2..#maxsub] do 
               S:=S meet maxsub[si];
             end for;
               J1 := VerticalJoin(
            Matrix(IntegerRing(), Morphism(S, Mp)), MatrixRing(Integers(),e)!p);
//Calculate the right idealiser of S
               R1 := BasisMatrix(Rowspace(J1));
               R1i:=Matrix(Rationals(),R1)^(-1);
               BB:=[]; 
               for i in [1..e] do
             Append(~BB,Matrix(Rationals(),R1)*Matrix(Rationals(),RM1[i])*R1i); 
               end for;
               SysEq:=KMatrixSpace(GF(p),e,e^2) ! 0;
               for  i in [1..e] do
               for  j in [1..e] do
               for  k in [1..e] do
                    SysEq[i][j+e*(k-1)] := GF(p) ! Integers() ! BB[i][j][k]; 
               end for;
               end for;
               end for;
               KK:=KernelMatrix(SysEq);
               if NumberOfRows(KK) gt 0 then 
                    //found submodule that has a bigger idealiser
                     J:=VerticalJoin(Matrix(IntegerRing(),KK),MatrixRing(Integers(),e) !p);
                     R:=BasisMatrix(RowSpace(J));
                     RS1:=(1/p*Matrix(Rationals(),R));
                     RS1i:=RS1^(-1);
                     RB1:=RS1*RB1;
//Calculate new regular representation
                   BB:=[];
                   for j in [1..e] do
                       Append(~BB, RS1 * Matrix(Rationals(),RM1[j]) * RS1i);
                   end for;

                   RM1:=[];
                   for j in [1..e] do
                       s:=MatrixRing(Rationals(),e) ! 0;
                       for k in [1..e] do s +:= RS1[j][k]*BB[k]; end for;
                       Append(~RM1,Matrix(Integers(),s));
                   end for;
                      new_order:=true;
               end if;
        end while;  // new order found
       "Number of Times through 1 :", number;
       "for Prime :", p;
     end for; //p in L

     det:=Determinant(RB1);
     "det,DFF*d^2"; det,DFF*det^2;
     DFFhereditary:=Integers() ! (DFF*det^2);
     "New DFF:", DFFhereditary;

     if Abs(DFFhereditary) eq 1 then
       return 1,realsi;
     end if;

     L := Factorization(DFFhereditary);
     "New Factorization:", L;

// we know that now the order is hereditary. 
// idealiser of a maximal ideal then yields an overorder. 
     for pt in L do
       p := pt[1];

       "Prime p:", p;
        number:=0;

        new_order := true;
        while new_order do
          number +:= 1;
        new_order:=false;
  LM1:=[]; //Left regular representation
       for i in [1..e] do X:=Parent(RM1[1]) ! 0;
                       for j,k  in [1..e] do
                       X[j][k]:=RM1[j][i][k];
               end for;
               Append(~LM1,X);
       end for;
           Ep := sub<MatrixRing(GF(p),e) | [ Matrix(GF(p),X) : X in RM1 cat LM1]>;
//take both right and left regular reps. since we are interested in twosided ideals
           Mp := RModule(Ep);
           maxsub := MaximalSubmodules(Mp);
           "#maxsub:", #maxsub;
            for si in [1..#maxsub] do 
               S:=maxsub[si];
            J1 := VerticalJoin(
            Matrix(IntegerRing(), Morphism(S, Mp)), MatrixRing(Integers(),e)!p);
//Calculate the right idealiser of S
               R1 := BasisMatrix(Rowspace(J1));
               R1i:=Matrix(Rationals(),R1)^(-1);
               BB:=[]; 
               for i in [1..e] do
             Append(~BB,Matrix(Rationals(),R1)*Matrix(Rationals(),RM1[i])*R1i); 
               end for;
               SysEq:=KMatrixSpace(GF(p),e,e^2) ! 0;
               for  i,j,k in [1..e] do
                    SysEq[i][j+e*(k-1)] := GF(p) ! Integers() ! BB[i][j][k]; 
               end for;
               KK:=KernelMatrix(SysEq);
               if NumberOfRows(KK) gt 0 then 
                    //found submodule that has a bigger idealiser
                     J:=VerticalJoin(Matrix(IntegerRing(),KK),MatrixRing(Integers(),e) !p);
                     R:=BasisMatrix(RowSpace(J));
                     RS1:=(1/p*Matrix(Rationals(),R));
                     RS1i:=RS1^(-1);
                     RB1:=RS1*RB1;
//Calculate new regular representation
                   BB:=[];
                   for j in [1..e] do
                       Append(~BB, RS1 * Matrix(Rationals(),RM1[j]) * RS1i);
                   end for;

                   RM1:=[];
                   for j in [1..e] do
                       s:=MatrixRing(Rationals(),e) ! 0;
                       for k in [1..e] do s +:= RS1[j][k]*BB[k]; end for;
                       Append(~RM1,Matrix(Integers(),s));
                   end for;
                      new_order:=true;
                      break si; 
               end if;
             end for; //si
        end while;  // new order found
       "Number of Times through 2 :", number;
       "for Prime :", p;
     end for; //p in L

     det:=Determinant(RB1);
     "det,DFF*d^2"; det,DFF*det^2;
     DFFmax:=Integers() ! (DFF*det^2);
     "New DFFmax:", DFFmax;

     if Abs(DFFmax) eq 1 then
       return 1,realsi;
     end if;

     L := Factorization(DFFmax);
     "New Factorization:", L;


SSP:=[];
     ss := 1; // it is the LCM of the local schurindices
     if realsi eq -1 then ss:=2; end if;
     for pt in L do
       p := pt[1];
       dp := Decomposition(K,p);
       normexponent := ExactQuotient(d, dp[1][2]); // d = deg(K)
       // So normexponent = number of primes over p times the exponent of the norm of one primideal over p 
       m2m := ExactQuotient(pt[2],  normexponent); // should be m(m-tp) 
       m2m:=ExactQuotient(m2m,rootofdimKofE);
       m2m:=-m2m + rootofdimKofE;
       ssp:=ExactQuotient(rootofdimKofE,m2m);
        Append(~SSP,[p,ssp]);
       ss:=LCM(ssp,ss);
     end for;
     return ss, SSP,realsi;
 end function;