/*
  Ulrich Thiel, 2020

  This Magma program computes Dirichlet units of an order using the proof of
  the Dirichlet unit theorem, made constructive in a brutal and inefficient way.

  The purpose is only to get a better understanding of the proof and to see that it is indeed constructive in principle.

  There may be some errors due to rounding (and other errors) but in principle
  it works, I think.

  Here's an example that (surprisingly) works:

  Attach("DirichletUnits.i.m");
  R<X> := PolynomialRing(Rationals());
  f := X^3 - X^2 + 4*X - 8;
  L := NumberField(f);
  O := MaximalOrder(L);
  time DirichletUnits(O);
    ...
    Time: 137.990

*/
intrinsic ArraySum(X) -> RngElt
{}
  if IsEmpty(X) then
    return 0;
  else
    return &+X;
  end if;
end intrinsic;

intrinsic DirichletUnits(O::RngOrd: delta:=0.000001)
{}

  n := #Basis(O);
  r := #RealEmbeddings(O.1);
  s := (#Conjugates(O.1) - r) div 2; // 1/2 the number of complex embeddings

  Sprintf("Embeddings: r=%o, s=%o", r, s);
  Sprintf("Rank of unit group is %o", r+s-1);

  //Conjugates has the real embeddings in front, but then pairs of
  //complex embeddings. I'm using numbering where one system of complex
  //embeddings comes next, then the conjugate system.
  //Function sigma yields this numbering.
  sigma := function(alpha, i)
    if i le r then
      return Conjugates(alpha)[i];
    elif i gt r and i le r+s then
      return Conjugates(alpha)[r + 2*(i-r) - 1];
    else
      return Conjugates(alpha)[r + 2*(i-r-s)];
    end if;
  end function;

  c := function(i)
    if i le r then
      return 1;
    else
      return 2;
    end if;
  end function;

  Lambda := MinkowskiLattice(O);

  //Inverse of Minkowski map j:O->Lambda
  jinverse := function(x)
    x:=Coordinates(Lambda,x);
    return &+[x[i]*O.i : i in [1..#x] ];
  end function;

  disc := Discriminant(O);
  print "Discriminant of O is "*Sprint(disc);

  pi := Pi(Parent(1.0));
  C := (2/pi)^s * Sqrt(Abs(disc)) + delta;
  Sprintf("Universal bound C is %.10o", C);

  eps := [];
  for i:=1 to r+s-1 do
    Sprintf("Constructing eps[%o]", i);
    IndentPush();
    k := 1;
    alphai := [One(O)];
    Sprintf("alpha[%o,1] = %o", i, alphai[1]);

    while true do
      k +:= 1;
      Ci := C;
      CiArr := [ Abs(sigma(alphai[k-1],j)) - delta : j in [1..r+s] ];
      CiArr[i] := Ci/&*[ CiArr[j]^c(j) : j in {1..r+s} diff {i}];
      Sprintf("C[%o,%o] = %.10o", i, i,Ci);
      CiAlpha := Sqrt(ArraySum([ CiArr[j]^2 : j in [1..r]]) + 4*ArraySum([ CiArr[j]^2 : j in [r+1..r+s]]));
      Sprintf("CiAlpha = %.10o", CiAlpha);
      Ei := ShortVectors(Lambda, CiAlpha^2 + delta);
      Sprintf("Lambda has %o vectors of norm <= CiAlpha", #Ei);
      Ei := [ jinverse(x[1]) : x in Ei ];
      Ei := [ x : x in Ei | Abs(Norm(x)) le Ci and { Abs(sigma(x, j)) lt Abs(sigma(alphai[k-1], j)) : j in {1..n} diff {i} } eq {true} and x notin alphai]; //last condition to avoid precision errors
      Sprintf("E[%o] contains %o vectors", i, #Ei);
      Append(~alphai, Ei[1]);
      Sprintf("alpha[%o,%o] = %o", i,k, alphai[k]);

      found := false;
      for j:=1 to k do
        for l:=j+1 to k do
          epsi := alphai[l]/alphai[j];
          epsiinv := alphai[j]/alphai[l];
          if epsi in O and epsiinv in O then
            Sprintf("Found eps[%o] = %o", i, O!epsi );
            found := true;
            break j;
          end if;
        end for;
      end for;

      if found then
        Append(~eps, epsi);
        break;
      end if;

    end while;

    IndentPop();

  end for;

  //Analyze cosets O^* / <eps>
  //eps := eps^2; //for testing
  print "Determining potential cosets:";
  IndentPush();
  jstar := function(alpha)
    return [ c(i)*Log(Abs(sigma(alpha,i))) : i in [1..n] ];
  end function;

  Cprime := 1/2 * &+[ Sqrt(&+[jstar(eps[i])[j]^2 : j in [1..n]]) : i in [1..r+s-1] ];

  C := Sqrt(n) * Exp(Cprime);

  Sprintf("C = %.10o", C);

  U := [ jinverse(x[1]) : x in ShortVectors(Lambda, C^2) | x[2] - n gt 10^-50 and Norm(jinverse(x[1])) in {1,-1}]; //the weird first condition should remove torsion units

  TU,f := TorsionUnitGroup(O); //yeah, I'm cheating, but we know how to compute
                              //torsion units in principle
  zeta := f(TU.1);
  o := Order(TU.1);
  for i:=1 to #eps do
    for k:=1 to o do
      Exclude(~U, zeta^k*O!eps[i]);
      Exclude(~U, zeta^k*(O!eps[i])^-1);
    end for;
  end for;


  print U;

end intrinsic;

intrinsic QuadraticFundamentalUnit(d::RngIntElt) -> SeqEnum
{}

  require d gt 0: "d must be positive";
  require IsSquarefree(d): "d must be squarefree";

  if d mod 4 in {2,3} then
    b := 0;
    while true do
      b +:= 1;
      f,a := IsSquare(b^2*d + 1);
      if f then
        return [a,b];
      end if;
      f,a := IsSquare(b^2*d - 1);
      if f then
        return [a,b];
      end if;
    end while;
  else
    b := 0;
    while true do
      b +:= 1;
      f,a := IsSquare(b^2*d + 4);
      if f then
        return [a,b];
      end if;
      f,a := IsSquare(b^2*d - 4);
      if f then
        return [a/2-b/2,b];
      end if;
    end while;
  end if;

end intrinsic;
