/* vim: set syntax=magma : */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// Stefano Marseglia, stefano.marseglia89@gmail.com                        //
// https://stmar89.github.io/index.html                                    //
//                                                                         //
// Code to compute the examples in the paper:                              //
// "Computing square-free polarized abelian varieties over finite fields"  //
// Mathematics of Computation 90 (2021), no. 328, 953–971.                 //
// https://doi.org/10.1090/mcom/3594                                       //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////

//////////////
// The packages are available at:
// https://github.com/stmar89/
//////////////

    AttachSpec("~/AlgEt/spec");
    AttachSpec("~/AlgEt/specMod");
    AttachSpec("~/AlgEt/specMtrx");
    AttachSpec("~/AbVarFq/spec");

//////////////////////////////////
//Example 7.2
//////////////////////////////////

    _<x>:=PolynomialRing(Integers());
    h:=x^4+2*x^3-7*x^2+22*x+121;
    AVh:=IsogenyClass(h);
    iso:=IsomorphismClasses(AVh);
    PHI:=pAdicPosCMType(AVh);

    for iA in [1..#iso] do
        A:=iso[iA];
        N:=0;
        repeat
            printf ".";
            N+:=1;
            test,pols_deg_N:=HasPolarizationsOfDegree(A,PHI,N);
        until test;
        for pol in pols_deg_N do
            aut:=#PolarizedAutomorphismGroup(pol);
        end for;
        printf "\nThe %o-th abelian variety of the list has %o polarizations of minimal degree %o, each with %o polarized automorphisms\n",iA,#pols_deg_N,N,aut;
    end for;

//////////////////////////////////
//Example 7.3
//////////////////////////////////

    _<x>:=PolynomialRing(Integers());
    h:=x^6-2*x^5-3*x^4+24*x^3-15*x^2-50*x+125;
    AVh:=IsogenyClass(h);
    "The isogeny class determined by the polynomial ",WeilPolynomial(AVh),"is defined over ",FiniteField(AVh)," and has dimension ",Dimension(AVh);
    S1:=ZFVOrder(AVh);
    O:=MaximalOrder(Algebra(S1));
    over_orders:=FindOverOrders(S1);
    for S in over_orders do
        "Is Gorenstein?", IsGorenstein(S),
        "\tIs conj. stable?", S eq ComplexConjugate(S),
        "\tIndex in O:",Index(O,S),
        "\tTors Units:",#TorsionSubgroup(UnitGroup(S));
    end for;
    
    iso:=IsomorphismClasses(AVh);
    "There are",#iso, "isomorphism classes.";
    
    for iA in [1..#iso] do
        A:=iso[iA];
        dir_sum:=DirectSumRepresentation(DeligneModule(A));
        I:=dir_sum[1,1];
        map:=dir_sum[1,2];
        assert map(1) eq 1;
        S:=EndomorphismRing(A);
        assert S eq MultiplicatorRing(I);
        "Is the ideal representing the",iA,"-th abelian variety invertible in its multiplicator ring?",IsInvertible(S!!I);
    end for;
        
    "We compute a CM-type with the prescribed p-adic positivity:";
    PHI:=pAdicPosCMType(AVh);
    
    for iA in [1..#iso] do
        A:=iso[iA];
        test,princ_pols:=IsPrincipallyPolarized(A,PHI);
        if test then
            aut:=#PolarizedAutomorphismGroup(princ_pols[1]);
            assert forall{ pol : pol in princ_pols | #PolarizedAutomorphismGroup(pol) eq aut};
        end if;
        "The",iA,"-th abelian variety has",#princ_pols,"principal polarizations, with",aut,"polarized automorphisms.";
    end for;

//////////////////////////////////
//Example 7.4
//////////////////////////////////

    _<x>:=PolynomialRing(Integers());
    h:=x^8-5*x^7+13*x^6-25*x^5+44*x^4-75*x^3+117*x^2-135*x+81;
    AVh:=IsogenyClass(h);
    "The isogeny class determined by the polynomial ",WeilPolynomial(AVh),"is defined over ",FiniteField(AVh)," and has dimension ",Dimension(AVh);
    S1:=ZFVOrder(AVh);
    O:=MaximalOrder(Algebra(S1));
    over_orders:=FindOverOrders(S1);
    for S in over_orders do
        "Is Gorenstein?", IsGorenstein(S),
        "\tIs conj. stable?", S eq ComplexConjugate(S),
        "\tIndex in O:",Index(O,S),
        "\tTors Units:",#TorsionSubgroup(UnitGroup(S));
    end for;
    
    "We compute the isomorphism classes with endomorphism ring which is conjugate stable";
    iso:=[ A : A in IsomorphismClasses(AVh) | ComplexConjugate(S) eq S where S:=EndomorphismRing(A)];
    "There are",#iso, "such isomorphism classes.";
    
           
    "We compute a CM-type with the prescribed p-adic positivity:";
    PHI:=pAdicPosCMType(AVh);
    
    for iA in [1..#iso] do
        A:=iso[iA];
        test,princ_pols:=IsPrincipallyPolarized(A,PHI);
        if test then
            aut:=#PolarizedAutomorphismGroup(princ_pols[1]);
            assert forall{ pol : pol in princ_pols | #PolarizedAutomorphismGroup(pol) eq aut};
        end if;
        "The",iA,"-th abelian variety has",#princ_pols,"principal polarizations, with",aut,"polarized automorphisms.";
    end for; 

    "we notice that there are abelian varieties with the same endomorphism ring whose groups of points are not isomorphic";
    iso:=IsomorphismClasses(AVh);
    for S in over_orders do if not IsGorenstein(S) then
        iso_S:=[A: A in iso|EndomorphismRing(A) eq S];
        group_of_points:={ ElementaryDivisors(RationalPoints(A)) : A in iso_S };
        if #group_of_points ne 1 then
            "here are the elementary divisors of the groups we find";
            group_of_points;
        end if;
    end if; end for;
    
//////////////////////////////////
//Example 8.1
//////////////////////////////////

    _<x>:=PolynomialRing(Integers());
    h:=(x^4-4*x^3+8*x^2-12*x+9)*(x^4-2*x^3+2*x^2-6*x+9);
    AVh:=IsogenyClass(h);
    "The isogeny class determined by the polynomial ",WeilPolynomial(AVh),"is defined over ",FiniteField(AVh)," and has dimension ",Dimension(AVh);

    "We compute the isomorphism classes with endomorphism ring which is conjugate stable";
    iso:=[ A : A in IsomorphismClasses(AVh) | ComplexConjugate(S) eq S where S:=EndomorphismRing(A)];
    "There are",#iso, "such isomorphism classes.";
    
    "We compute a CM-type with the prescribed p-adic positivity:";
    PHI:=pAdicPosCMType(AVh);
    
    for iA in [1..#iso] do
        A:=iso[iA];
        test,princ_pols:=IsPrincipallyPolarized(A,PHI);
        if test then 
            for pol in princ_pols do
                PeriodMatrix(pol,PHI);
            end for;
        end if;
    end for;