/* vim: set syntax=magma :*/

/////////////////////////////////////////////////////
// Stefano Marseglia, stefano.marseglia89@gmail.com
// https://stmar89.github.io/index.html
/////////////////////////////////////////////////////

/*
    Here you can find the code to reproduce the examples in the paper:
        "Modules over orders, conjugacy classes of integral matrices, and abelian varieties over finite fields",
        Stefano Marseglia
        https://arxiv.org/abs/2208.05409
    One needs several functionalities for étale algebra.
    Please clone the following repository on your machine:
        https://github.com/stmar89/AlgEt

    WARNING:
    The isomorphism testing (IsIsomorphic) and the computation of the isomorphism classes (IsomorphismClasses) 
    with the Method Magma is based on results from 
        "The conjugacy problem in GL(n, Z)", T. Hofmann, B. Eick and E. A. O’Brien, J. Lond. Math. Soc. (2) 100 (2019), no. 3, 731–756.
    and the corresponding implementation of their code, which is included in Magma.


    In the paper, we suggest to use instead to use the results from
        "Computation of lattice isomorphisms and the integral matrix similarity problem", T. Hofmann, W. Bley and H. Johnston, Forum Math. Sigma, 10:Paper No. e87, 36, 2022
        https://arxiv.org/abs/2202.03526
    Their code at the moment is not implemented in MAGMA, but only in the julia package Hecke/Nemo. 
    Requires julia 1.7.0 or higher and Hecke http://www.thofma.com/Hecke.jl/dev/ 

    Here are the instructions to use the hybrid Magma+julia option (on Unix):
    - install julia
    - add julia to the PATH variable
        eg. modify and add the line "PATH="$PATH:/home/username/julia-1.7.3/bin/julia-1.7.3/bin"; export PATH" to your .bashrc
            then run "source .bashrc" or restart your system.
            to make sure: the following command in Magma, should not trigger an error
            > Pipe("julia","");
    - open julia and run the following code (the last command might take a while, but it will make every future startup of julia+Hecke much faster).
    julia> using Pkg
    julia> Pkg.add("Hecke")
    julia> using Hecke
    julia> Hecke.build()
    - Take note of the path which is printed out (something like "/tmp/Hecke.so" or "/tmp/Hecke.dylib") and change it in the variable "julia_compiled" below.
    
    The code was last tested on May 30, 2024, on Magma V2.28-8, julia 1.7.3, Hecke 0.14.13.

*/

// Before starting with the examples, modify the following line and attach the packages AlgEt and AlgEtMod

path:="~/AlgEt/"; // Put here the path where you have cloned the repo AlgEt
julia_compiled:="/tmp/Hecke.so"; //modify this according to the output of Hecke.build, as explained above.
AttachSpec(path cat "spec");
AttachSpec(path cat "specMod");

// ###################
// ### Example 6.1 ###
// ###################

_<x>:=PolynomialRing(Integers());
m1:=x^2-x+3;
m2:=x^2+x+3;
s1:=2;
s2:=1;
h:=m1^s1*m2^s2; h;
q:=Integers() ! Truncate( ConstantCoefficient(h)^(2/Degree(h)) );

K1:=NumberField(m1);
K2:=NumberField(m2);
K:=EtaleAlgebra([K1,K2]); // K = K1 x K2
pi:=PrimitiveElement(K);
R:=Order([pi,q/pi]);
O:=MaximalOrder(K);
V:=EtaleAlgebra([K1,K1,K2]); // V = K1^s1 x K2^s2
m:=NaturalAction(K,V); // m:K -> V component-wise diagonal action of K on V

// To run the Magma-only version of the code, use the following line:
time classes:=IsomorphismClasses(R,m); //~1 minute

// To use the hybrid Magma-julia version, follow the instrucions above and then run
time classes:=IsomorphismClasses(R,m : Method:="julia -J " cat julia_compiled cat " " cat path cat "AlgEtQMod/");
                                                                                                   // It should take ~12 seconds.
#classes;
// classes now contains the representatives of the isomorphism classes of the sub-R-modules of V.
// By Theorem 5.1, they correspond to the Fq-isomorphism classes of abelian varieties in the isogeny class determined by h.



// ###################
// ### Example 6.2 ###
// ###################
// We want to find the representatives of the conjugacy classes of integral matrices with minimal polynomial m1*m2 and characteristic polynomial h.
// This code is essentially the same of the one from the previous example, but replacing R=Z[pi,q/pi] with R=Z[pi]. Since E is smaller than R, there will be more candidates. So this computation will take more time to run than the previous one.

_<x>:=PolynomialRing(Integers());
m1:=x^2-x+3;
m2:=x^2+x+3;
s1:=2;
s2:=1;
h:=m1^s1*m2^s2; h;
K1:=NumberField(m1);
K2:=NumberField(m2);
K:=EtaleAlgebra([K1,K2]); // K = K1 x K2
pi:=PrimitiveElement(K); 
R:=Order([pi]);
O:=MaximalOrder(K);
V:=EtaleAlgebra([K1,K1,K2]); // V = K1^s1 x K2^s2
m:=NaturalAction(K,V); // m:K -> V component-wise diagonal action of K on V


// To run the Magma-only version of the code, use the following line:
time classes:=IsomorphismClasses(R,m); // Attention! This is very slow!

// To use the hybrid Magma-julia version, follow the instrucions above and then run
time classes:=IsomorphismClasses(R,m : Method:="julia -J " cat julia_compiled cat " " cat path cat "AlgEtQMod/");
                                                                                                   // It should take ~200 secs
#classes;
// By Theorem 4.1, they correspond to the conjugacy classes of matrices with with minimal polynomial m1*m2 and characteristic polynomial h.
// We now print the 4 matrices. These matrices might not be the same one as in the paper, but they are conjugate to them.
for I in classes do
   mat:=Matrix(AbsoluteCoordinates([m(pi)*z : z in ZBasis(I)],ZBasis(I)));
   assert MinimalPolynomial(mat) eq m1*m2;
   assert CharacteristicPolynomial(mat) eq h;
   printf "%o\n\n",mat;
end for;

quit;