// Examples for the package 'ellfib'.
// GDB, DJR 19 June 2008

// Lines of input code below - expected output follows the input
// but commented out (either by // ... or by /* ... */).

// cd to the directory containing the file 'ellfib.mag'.
// Run magma.
// Paste in the following lines:

Attach("ellfib.mag");	// this loads the ellfib package

// (5.1) An Halphen fibration with mu=2
k := Rationals();
P3<x,y,z,t> := ProjectiveSpace(k,3);
X := Scheme(P3,t^3 - x^3 + y^2*z + 2*x*z^2 - z^3);
IsNonsingular(X);
  // true
O := X ! [0,1,0,0];
R := X ! [1,0,1,0];
L := Scheme(P3,[y,t]);
PandR := Intersection(X,L);
P := [ Z : Z in IrreducibleComponents(PandR) | Degree(Z) eq 2 ][1];  P;
  // Scheme over Rational Field defined by x^2 + x*z - z^2, y, t
A2 := LinearSystem(P3,2);
H0 := ImposeBasepoint(X,A2,P,2);
H := ImposeBasepoint(X,H0,O,2);
H;
  // Linear system on Projective Space of dimension 3
  // with 2 sections: x^2 + x*z - z^2, t^2
C := Curve(Intersection(X, Scheme(P3, t^2 + x^2 + x*z - z^2)));
assert IsIrreducible(C);
Genus(C);
  // 1

// (5.2) Geiser and Bertini involutions
X := Scheme(P3,x^3 + y^3 + z^3 + 3*t^3);
P := X ! [1,1,1,-1];
iP := GeiserInvolution(X,P);
DefiningEquations(iP);
  /*
  [
      -x*y + y^2 - x*z + z^2 - 3*x*t - 3*t^2,
      x^2 - x*y - y*z + z^2 - 3*y*t - 3*t^2,
      x^2 + y^2 - x*z - y*z - 3*z*t - 3*t^2,
      -x^2 - y^2 - z^2 - x*t - y*t - z*t
  ]
  */
TP := TangentSpace(X,P);
CP := Curve(Intersection(X,TP));
iP(CP);
  // Scheme over Rational Field defined by z + t, y + t, x + t
Support(iP(CP));
  // { (-1 : -1 : -1 : 1) }
L := Scheme(P3,[x-y,z+t]);
XL := Intersection(X,L);
Q := [ Z : Z in IrreducibleComponents(XL) | Degree(Z) eq 2 ][1];
iQ := BertiniInvolution(X,Q);
DefiningEquations(iQ);
  /*
  [
    x^2*y^3 - 5/6*x*y^4 + 5/6*y^5 - 1/6*x^2*y^2*z - 1/6*x*y^3*z - 2/3*x^2*y*z^2 
        - 2/3*y^3*z^2 + x^2*z^3 - 2/3*x*y*z^3 + 11/6*y^2*z^3 - 1/3*x*z^4 + y*z^4
        - 2/3*z^5 - 1/2*x^2*y^2*t - 1/2*x*y^3*t + 3*x^2*z^2*t - 1/2*x*y*z^2*t + 
        9/2*y^2*z^2*t - 2*x*z^3*t + 3/2*y*z^3*t - 2*z^4*t + 2*x^2*y*t^2 + 
        2*y^3*t^2 + 3*x^2*z*t^2 - 3/2*x*y*z*t^2 + 7/2*y^2*z*t^2 - 4*x*z^2*t^2 - 
        3*y*z^2*t^2 + x^2*t^3 - 10/3*x*y*t^3 + 5/6*y^2*t^3 - 10/3*x*z*t^3 - 
        5/2*y*z*t^3 + 16/3*z^2*t^3 - x*t^4 + 3*y*t^4 + 6*z*t^4 + 2*t^5,
    -5/6*x^2*y^3 + 5/6*x*y^4 - y^5 - 1/6*x^2*y^2*z + 1/6*y^4*z - 2/3*x*y^2*z^2 +
        2/3*y^3*z^2 + x^2*z^3 + 1/6*x*y*z^3 + x*z^4 - 1/6*y*z^4 - 1/2*x^2*y^2*t 
        + 1/2*y^4*t + 9/2*x^2*z^2*t - 1/2*x*y*z^2*t + 3*y^2*z^2*t + 3/2*x*z^3*t 
        - 3/2*y*z^3*t - 2*z^4*t + 2*x*y^2*t^2 - 2*y^3*t^2 + 7/2*x^2*z*t^2 - 
        3/2*x*y*z*t^2 + 3*y^2*z*t^2 - 3*x*z^2*t^2 - 4*y*z^2*t^2 - 2*z^3*t^2 - 
        5/3*x^2*t^3 - 5/6*x*y*t^3 - 2*y^2*t^3 - 5/2*x*z*t^3 - 17/6*y*z*t^3 + 
        22/3*z^2*t^3 + 3*x*t^4 + 1/2*y*t^4 + 6*z*t^4 - 4*t^5,
    -2*x^2*y^2*z + 5/3*x*y^3*z - 5/3*y^4*z + 1/3*x^2*y*z^2 + 1/3*x*y^2*z^2 + 
        2/3*x^2*z^3 + 2/3*y^2*z^3 + 1/3*x*z^4 - 4/3*y*z^4 - z^5 - 3*x^2*y^2*t + 
        5/2*x*y^3*t - 5/2*y^4*t + 3/2*x^2*y*z*t + 3/2*x*y^2*z*t - x*z^3*t - 
        7/2*y*z^3*t - 4*z^4*t + 3/2*x^2*y*t^2 + 3/2*x*y^2*t^2 - 6*x^2*z*t^2 - 
        6*y^2*z*t^2 - 3/2*x*z^2*t^2 - 3/2*y*z^2*t^2 - 5/2*z^3*t^2 - 20/3*x^2*t^3
        - 20/3*y^2*t^3 + 43/6*x*z*t^3 + 13/6*y*z*t^3 + 9*x*t^4 + 3/2*y*t^4 - 
        23/2*z*t^4 - 13*t^5,
    x^2*y^2*z - 5/6*x*y^3*z + 5/6*y^4*z - 1/6*x^2*y*z^2 - 1/6*x*y^2*z^2 - 
        1/3*x*z^4 + 1/2*y*z^4 - 1/3*z^5 + 2*x^2*y^2*t - 5/3*x*y^3*t + 5/3*y^4*t 
        - 5/6*x^2*y*z*t - 5/6*x*y^2*z*t + 2/3*x^2*z^2*t + 2/3*y^2*z^2*t - 
        1/2*x*z^3*t + 7/6*y*z^3*t - 5/6*z^4*t - x^2*y*t^2 - x*y^2*t^2 + 
        4*x^2*z*t^2 + 4*y^2*z*t^2 - 1/2*x*z^2*t^2 - 1/2*y*z^2*t^2 - 2*z^3*t^2 + 
        14/3*x^2*t^3 + 14/3*y^2*t^3 - 5*x*z*t^3 - 5/2*y*z*t^3 - 11/6*z^2*t^3 - 
        19/3*x*t^4 - 4/3*y*t^4 + 20/3*z*t^4 + 9*t^5
  ]
  */

// (5.3) Eckardt points
X := Scheme(P3, x^3 + y^3 + z^3 + 2*t^3);
P := X ! [1,-1,0,0];
iP := GeiserInvolution(X,P);
[ f div GCD(E) : f in E ] where E is DefiningEquations(iP);
  // [ y, x, z, t ]

// (5.4) An example of untwisting

f1 :=
57645*x^2*y^3 + 47234*x*y^4 - 9963*y^5 + 23490*x^2*y^2*z + 97322*x*y^3*z +
70056*y^4*z - 26730*x^2*y*z^2 - 33603*x*y^2*z^2 + 5751*y^3*z^2 +
47925*x^2*z^3 + 85664*x*y*z^3 - 5373*y^2*z^3 + 41480*x*z^4 + 72990*y*z^4 +
4095*z^5 + 8100*x^2*y^2*t + 157516*x*y^3*t + 148392*y^4*t - 200880*x^2*y*z*t -
25896*x*y^2*z*t + 182664*y^3*z*t + 9720*x^2*z^2*t - 10800*x*y*z^2*t -
42408*y^2*z^2*t + 118912*x*z^3*t + 194220*y*z^3*t + 109800*z^4*t -
124740*x^2*y*t^2 - 27990*x*y^2*t^2 + 96462*y^3*t^2 - 42120*x^2*z*t^2 -
112938*x*y*z*t^2 - 70722*y^2*z*t^2 + 24042*x*z^2*t^2 + 28314*y*z^2*t^2 +
63558*z^3*t^2 + 118530*x^2*t^3 + 111736*x*y*t^3 - 48186*y^2*t^3 +
157684*x*z*t^3 + 176616*y*z*t^3 + 14958*z^2*t^3 + 247316*x*t^4 + 338796*y*t^4
+ 265536*z*t^4 + 123444*t^5;
f2 :=
20232*x^2*y^3 + 27216*x*y^4 + 6600*y^5 - 66429*x^2*y^2*z - 29187*x*y^3*z +
40250*y^4*z + 25596*x^2*y*z^2 - 8532*x*y^2*z^2 - 42800*y^3*z^2 +
24507*x^2*z^3 + 23436*x*y*z^3 + 3585*y^2*z^3 - 4185*x*z^4 + 35420*y*z^4 -
38240*z^5 - 48978*x^2*y^2*t + 77706*x*y^3*t + 128092*y^4*t - 84456*x^2*y*z*t -
85428*x*y^2*z*t - 11724*y^3*z*t + 65322*x^2*z^2*t + 26676*x*y*z^2*t -
8214*y^2*z^2*t + 100710*x*z^3*t + 125152*y*z^3*t + 25500*z^4*t -
196596*x^2*y*t^2 - 75438*x*y^2*t^2 + 122086*y^3*t^2 - 106596*x^2*z*t^2 -
104598*x*y*z*t^2 + 366*y^2*z*t^2 + 4590*x*z^2*t^2 - 6786*y*z^2*t^2 +
144574*z^3*t^2 - 62424*x^2*t^3 - 63612*x*y*t^3 - 16932*y^2*t^3 -
105030*x*z*t^3 + 1972*y*z*t^3 - 98056*z^2*t^3 + 117720*x*t^4 + 231884*y*t^4 +
36888*z*t^4 + 247412*t^5;
P1 := ProjectiveSpace(k,1);
f := map< P3 -> P1 | [f1,f2] >;
  // Note: since involutions are returned as maps P3 - -> P3, we need
  // to define f as a map P3 - -> P1 in order to check the untwisting.
time existence, Q := HasMaximalCentre(f,X); assert existence;
  // Time: 64.190
Q;
  // Scheme over Rational Field defined by
  // z^2 - 31/4*z*t - 5/4*t^2,  x + 3/2*z + 3/2*t,  y - 3/2*z - 1/2*t
Degree(Q);
  // 2
k2<w> := Degree2SplittingField(Q);
Support(Q,k2);
  // { (w : -w - 1 : 1/3*(-2*w - 3) : 1),
  // (1/8*(-8*w - 117) : 1/8*(8*w + 109) : 1/12*(8*w + 105) : 1) }
iQ := BertiniInvolution(X,Q);
g := iQ * f;
Eg := DefiningEquations(g);
assert IsDivisibleBy(x*Eg[2] - y*Eg[1], DefiningEquation(X));

// (5.5) The problem of minimality
X := Scheme(P3,x*t^2 + x^2*y + y^3 - z^3);
Z := Scheme(P3,[x,t,y^2+y*z+z^2]);
L1 := ImposeBasepoint(X, LinearSystem(P3,5), Z, 6);
L2 := Complement(L1,X);
#Sections(L2);
  // 5
P4<[a]> := ProjectiveSpace(k,4);
f := map< P3 -> P4 | Sections(L2) >;
f(X);
  /*
  Scheme over Rational Field defined by
  a[2]*a[4] - a[3]*a[4] + a[1]*a[5],
  a[1]*a[2]*a[4] - a[4]^3 - a[2]^2*a[5] - a[2]*a[3]*a[5] - a[3]^2*a[5],
  a[1]^2*a[2] + a[2]^3 - a[3]^3 - a[1]*a[4]^2
  */