read("things.gp");
read("h1.gp");
\\ This function returns the intersection number of the lines
\\ L^a_{mn} and L^b_{ij} where a and b are values 0-2 representing
\\ (123), (231), (312) respectively, and m,n,i,j are values in 1,3,5,7.
intersect(a,m,n,b,i,j) =
{
if(a == b,
if((m == i) && (n == j),
-2,
if((m == i) || (n == j), 1, 0)
),
if(Mod(b,3) == Mod(a+1,3),
if(Mod(m-n,8) == Mod(i+j,8),1,0),
if(Mod(i-j,8) == Mod(m+n,8),1,0)
)
)
}
\\ We define the matrix Q to be the matrix of intersection numbers of
\\ the 48 lines. The basis is ordered thus:
\\ L^123_11, L^123_13, ... , L^123_57, L^123_77, L^231_11, ...
Q = matrix(48,48);
{
for(a=0,2,
for(m=0,3,
for(n=0,3,
for(b=0,2,
for(i=0,3,
for(j=0,3,
Q[1+a*16+m*4+n,1+b*16+i*4+j] =
intersect(a,2*m+1,2*n+1,b,2*i+1,2*j+1)
)
)
)
)
)
)
}
\\ A0 = Z-basis for kernel of Q, that is, those combinations of lines
\\ which are principal divisors. This is because the kernel is Q consists
\\ of those sums of lines numerically equivalent to 0 (by definition).
\\ But on a K3 surface, numerical equivalence is the same as linear
\\ equivalence.
A0 = matrixqz(matker(Q),-2);
\\ P = plane section ( = Z-basis for positive definite part)
P = [1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]~;
\\ B = Q-basis for negative definite part
\\ so now the columns of A0, B and P together form a basis for the
\\ vector space spanned by the 48 lines.
B = matintersect(vecextract(matsupplement(A0),"29..48"),matker(Mat(P~*Q)));
\\ C = matrix which projects onto the subspace we want. We project
\\ onto the space spanned by P and B, such that the space spanned by
\\ the columns of A0 goes to 0. At the moment this is only about vector
\\ spaces.
C1 = concat(P,concat(B,A0));
C2 = matdiagonal([1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]);
C = C1*C2*C1^-1;
\\ Find a Z-basis for the lattice generated by the columns of C; this is
\\ a basis for Pic V. Why? Because the columns of C are the images of
\\ the 48 lines under C.
D = mathnf(C);
\\ E = matrix which maps NS(V) onto Z^20 (a left inverse for D)
E = matleftinverse(D) * C;
\\ Now we can use E and D to map between the space generated by the
\\ 48 lines and a 20-dimensional space representing Pic V.
\\ ---------------------------------------------------------------------
\\ From now on, we define functions which often take a group H. This
\\ represents a subgroup of our group Gamma, which is
\\ Gal(Q(e,t_1,t_2,t_3)/Q(t_1^4,t_2^4,t_3^4))
\\ where t_i are transcendental elements, and e is a primitive eighth
\\ root of unity. We fix a set of generators for Gamma, which in order
\\ are:
\\ g_1: t_1 -> i t_1
\\ g_2: t_2 -> i t_2
\\ g_3: t_3 -> i t_3
\\ g_4: \sqrt{2} -> -\sqrt{2}
\\ g_5: i -> -i
\\ where each generator fixes those elements not mentioned in its
\\ definition. The first four generate an Abelian group of type
\\ [4,4,4,2], which we call Gamma'. The last one doesn't commute with
\\ the others: in fact conjugating by the last generator inverts all
\\ the other, as in a dihedral group.
\\ The PARI object H is a vector with two elements. The first element
\\ is a 4x4 matrix. The columns of this matrix represent generators for
\\ the intersection of H with Gamma', expressed in terms of the first four
\\ generators given above. The second element of H is a vector. If H
\\ is contained in Gamma', it is an empty vector. Otherwise, it defines
\\ an element which, multiplied by the fifth generator of G given above,
\\ gives a fifth generator of H.
\\ For the remaining part of this file, we will use a running example to
\\ illustrate the algorithms. This example corresponds to the surface
\\ described in Section 5.3 of the thesis ("Example with a fibration").
\\ That surface is defined by
\\ X_0^4 + X_1^4 = 6X_2^4 + 12X_3^4
\\ so its coefficients are (1,1,-6,-12). The corresponding subgroup H
\\ is represented by the matrix-vector pair
\\ [4 0 0 0] [0]
\\ [ [0 2 1 1] , [1] ] written out nicely here for convenience.
\\ [0 0 1 0] [1]
\\ [0 0 0 1] [0]
\\ To see how to get these values from the coefficients, look at the
\\ function find_group below.
\\ Thus the group H is generated by the following five elements:
\\ g_1^4 (which is just the identity)
\\ g_2^2 (taking t_2 -> -t_2)
\\ g_2.g_3 (taking t_2 -> it_2 and t_3 -> it_3)
\\ g_2.g_4 (taking t_2 -> it_2 and \sqrt{2} -> -\sqrt{2}
\\ g_5.g_2.g_3 (taking t_2 -> -it_2 and t_3 -> -it_3 and i -> -i)
\\ ---------------------------------------------------------------------
\\ Return a code for the constant field fixed by H:
\\ 1 Q
\\ 2 Q(i)
\\ 3 Q(\sqrt(2))
\\ 4 Q(\sqrt(-2))
\\ 5 Q(i,\sqrt(2))
\\ In our example, H fixes both \sqrt{2} and i, so the constant field
\\ is Q.
constant_field(H) =
{
if(H[2] == [], \\ is i fixed?
if((H[1][4,] % 2) == [0,0,0,0], 5, 2), \\ yes - either Q(i) or Q(e)
if((H[1][4,] % 2) == [0,0,0,0], \\ no - look at the action
if((H[2][4] % 2) == 0, 3, 4), \\ on \sqrt{2}
1 \\ Both fixed - must be Q.
)
)
}
\\ ---------------------------------------------------------------------
\\ Define a couple of vectors because we don't want to have to write
\\ them out again. X4[i][j] is I^{(i-1)(j-1)}. Remember that indices
\\ in PARI start at 1. X2[i][j] is (-1)^{(i-1)(j-1)}.
X4 = [[1,1,1,1],[1,I,-1,-I],[1,-1,1,-1],[1,-I,-1,I],[1,1,1,1]];
X2 = [[1,1],[1,-1],[1,1]];
\\ Now we're concerned with finding the fixed field of a given subgroup,
\\ as described in Algorithm 3.16 of the thesis.
\\ Let A = Q[t_1^4, t_2^4, t_3^4] and B = Q(e)[t_1, t_2, t_3]. Then
\\ B is a Galois algebra over A, with group Gamma. We have a subgroup
\\ H of Gamma, and want to find the part of B fixed by H. This depends
\\ on the action of H on the elements
\\ x = (t_1)^{m_1} (t_2)^{m_2} (t_3)^{m_3} \sqrt{2}^{m_4}
\\ of which there are 128, independent over A, which generate B as
\\ an A[Gamma]-module. Each such x is taken to one of x, -x, ix, -ix
\\ by any given element of Gamma.
\\ Return a 128-element vector which contains the effect of applying
\\ the group element g to each of the 128 basis vectors. Each element
\\ of the result is one of 1, -1, i, -i.
doelt(g) =
{
local(A,B,C,D);
A = X4[g[1]+1];
B = X4[g[2]+1];
C = X4[g[3]+1];
D = X2[g[4]+1];
tensor(A, tensor(B, tensor(C, D)))[1,];
}
\\ Return a matrix which gives a set of generators for the part of B fixed
\\ by H, as an algebra over A. The matrix represents such a generating set
\\ on the (multiplicative) basis ((1+i), \sqrt(2), t_1, t_2, t_3). See
\\ Lemmas 3.13--3.15.
fixed_field(H) =
{
local(M, n, b, c, N);
\\ First find a basis for the fixed part as a module over A
n = matsize(H[1])[2];
M = matrix(n, 128);
\\ Do the generators which fix i first
for(i = 1, n, M[i,] = doelt(H[1][,i]));
\\ For our running example, we look at the action of each generator of H'
\\ on the elements x (as defined above).
\\ g_1^4 fixes everything (because it's the identity);
\\ g_2^2 fixes x if and only if m_2 is even;
\\ g_2.g_3 fixes x if and only if m_2+m_3 = 0 (mod 4)
\\ g_2.g_4 fixes x if and only if 2m_2 + m_4 = 0 (mod 4)
\\ Find out which basis elements were fixed by all of these
b = vector(128);
for(j = 1, 128,
b[j] = 1;
for(i = 1, n, if(M[i,j] != 1, b[j] = 0; break));
);
\\ Now b contains entries either 1 or 0, indicating which basis
\\ elements are fixed by H'.
\\ In our running example, we see that the elements fixed by all
\\ four generators of H' are those with (m_2=m_3=m_4=0) or (m_2=m_3=1
\\ and m_4=1). That is, x = t_1^{m_1} or x = \sqrt{2}t_1^{m_1}t_2^2t_3^2.
\\ These elements generate the fixed part of B as a module over A[i].
\\ See what happens when we apply the generator which takes i to -i
if(H[2] != [],
c = doelt(H[2]);
for(j = 1, 128,
if(b[j] == 1,
if(c[j] == I, b[j] = 1-I,
if(c[j] == -1, b[j] = I,
if(c[j] == -I, b[j] = 1+I)
)
)
)
)
);
\\ The vector b now indicates (by Lemma 3.14) a set of generators
\\ for the fixed part of B as a module over A.
\\ In our example, we are looking at the element g_5.g_2.g_3. It
\\ fixes both t_1^{m_1} and \sqrt{2}t_1^{m_1}t_2^2t_3^2, so these
\\ two families of elements generate the fixed part of B as a module
\\ over A.
\\ Next express our basis elements multiplicatively in terms of (1+i),
\\ \sqrt(2), t_1, t_2, t_3
N = matrix(5,129);
for(i = 1, 128,
if(b[i] != 0,
N[3,i] = (i - 1) % 4;
N[4,i] = ((i - 1) \ 4) % 4;
N[5,i] = ((i - 1) \ 16) % 4;
N[2,i] = (i - 1) \ 64;
if(b[i] == 1+I, N[1,i] = 1,
if(b[i] == I, N[1,i] = 2,
if(b[i] == 1-I, N[1,i] = 3)
)
)
)
);
\\ If i is fixed, stick that in too
if(H[2] == [], N[1,129] = 2);
\\ For our example, the matrix N has eight non-zero columns, corresponding
\\ to the eight generators of the fixed part of B. The columns for
\\ t_1^{m-1} look like [0, 0, m_1, 0, 0]~ and those for
\\ \sqrt{2}t_1^{m_1}t_2^2t_3^2 look like [0, 1, m_2, 2, 2]~.
\\ Now, we only want generators for the fixed part of B as an algebra
\\ over A, rather than as a module. It's clear that we don't need
\\ all four elements t_1^{m_1} for m_1=0,1,2,3; the single element
\\ t_1 will generate all of these. Similarly for the other four
\\ generators. By putting N into Hermite normal form, we find a minimal
\\ generating set.
\\ Find a set of multiplicative generators
N = mathnf(concat(N, matdiagonal([4,2,4,4,4])));
\\ Now N has two non-trivial columns: [0, 0, 1, 0, 0]~ and
\\ [0, 1, 0, 2, 2]~. These correspond to t_1 and \sqrt{2}t_2^2t_3^2,
\\ which do indeed generate the fixed part of B as an algebra over A.
\\ The whole of N is:
\\ [4 0 0 0 0]
\\ [0 2 0 0 1]
\\ N = [0 0 1 0 0]
\\ [0 0 0 4 2]
\\ [0 0 0 0 2]
}
\\ This is a purely cosmetic function that takes the matrix from
\\ the preceding function and turns it into a vector of monomials
\\ representing generators.
\\ For our example, it returns the PARI vector [t1, t3^2*t2^2*r2].
ff2gens(N) =
{
local(n, u, v, w, m);
n = matsize(N)[2];
v = vector(n);
u = [1 + I, 'r2, 't1, 't2, 't3];
w = [4, 2, 4, 4, 4];
for(i = 1, n,
v[i] = 1;
for(j = 1, 5,
v[i] *= u[j]^(N[j,i] % w[j])
);
if(N[1,i] == 2 || N[1,i] == 3, v[i] /= 2);
);
m = 0;
for(i = 1, n, if(v[i] != 1, m++));
w = vector(m);
m = 0;
for(i = 1, n, if(v[i] != 1, m++; w[m] = v[i]));
w;
}
\\ ---------------------------------------------------------------------
\\ Using the generators for the fixed part of B, we may easily
\\ test whether a few things are squares or not.
\\ Return a value:
\\ 1 a0a1a2a3 = square
\\ 2 a0a1a2a3 = 2 * square
\\ 3 a0a1a2a3 = -square
\\ 4 a0a1a2a3 = -2 * square
\\ 0 otherwise
\\ In our example, we find that \sqrt{2}(t_1t_2t_3)^2 is fixed by H;
\\ this means that for this subgroup to come from a surface, we must
\\ have 2a_0a_1a_2a_3 square. The function returns 2.
square_type(H) =
{
local(N);
N = fixed_field(H);
if(matsolvemod(N, [2,1,4,4,4], [0,0,2,2,2]~) == 0, 0,
if(matsolvemod(N, [4,2,4,4,4], [0,0,2,2,2]~), 1,
if(matsolvemod(N, [4,2,4,4,4], [0,1,2,2,2]~), 2,
if(matsolvemod(N, [4,2,4,4,4], [2,0,2,2,2]~), 3,
if(matsolvemod(N, [4,2,4,4,4], [2,1,2,2,2]~), 4, -1)
)
)
)
);
}
\\ ---------------------------------------------------------------------
\\ A little function to take out powers of 4. Only sensible for
\\ small numbers - but that's all we're using.
remove4s(n) = local(f); f = factor(n); f[,2] %= 4; factorback(f);
\\ This function implements Algorithm 3.17, which takes a subgroup H
\\ of Gamma and gives conditions on the coefficients a_i for the surface
\\ defined by them to be associated to a subgroup of H.
\\ All we're doing here is rearranging a few equations. We assume that
\\ the base field is Q. Each column of N will represent an equation:
\\ specifically, the column [a,b,c,d,e,f,g,h]~ means
\\ c_1^a c_2^b c_3^c = (1+i)^d \srqt{2}^e t_1^f t_2^g t_3^h
\\ and each of these is rearranged to give an expression for t_i in terms
\\ of everything else. The fact that fixed_field returns a matrix in
\\ Hermite normal form makes sure that everything works.
\\ In our running example, we have seen that t_1 and \sqrt{2}t_2^2t_3^2
\\ are fixed by the group H. We are looking for conditions on the
\\ coefficients [a_0, a_1, a_2, a_3] such that the Galois group
\\ Gal(Q(e, q_1, q_2, q_3)/Q)
\\ is contained in H, where q_i is a fourth root of (a_i/a_0) and we
\\ associate that Galois group to Gamma by associating t_i to q_i.
\\ The conditions will be (looking at the last three columns of N above)
\\ q_1 is rational
\\ q_2^4 is rational (a vacuous condition)
\\ \sqrt{2} q_2^2 q_3^2 is rational
\\ To solve this, we take a_0=1 and set each rational expression equal to a
\\ constant:
\\ q_1 = c_1
\\ q_2^4 = c_2
\\ \sqrt{2} q_2^2 q_3^2 = c_3
\\ and then solve for the a_i by taking appropriate powers:
\\ a_1 = c_1^4
\\ a_2 = c_2
\\ 2a_2a_3 = c_3^2
\\ Rearranging and modulo fourth powers, we get
\\ [a_0, a_1, a_2, a_3] = [1, 1, c_2, 2 c_3^2 c_2^3] .
example(H) =
{
local(N, v, n, c);
\\ First form the matrix N
N = concat([0,0,0;0,0,0;1,0,0;0,1,0;0,0,1],fixed_field(H)~)~;
\\ Rearrange the equations
for(i=3,5,
N[,i] *= (4 / N[i+3,i]);
for(j=3, i-1, N[,i] -= (N[j+3,i]/N[j+3,j])*N[,j]);
);
\\ Now form a monomial from each column
v = vector(4);
v[1] = 1;
for(i=3, 5,
v[i-1] = 'c1^(N[1,i] % 4) * 'c2^(N[2,i] % 4) * 'c3^(N[3,i] % 4);
n = (1+I)^N[4,i] * ('r2^N[5,i] % ('r2^2-2));
\\ Remove any unnecessary constants
c = gcd(real(content(n)), imag(content(n)));
n /= (c / remove4s(c));
v[i-1] *= n;
);
v;
}
\\ ---------------------------------------------------------------------
\\ Now we move onto the action on the Picard group. As pointed out in
\\ Section 3.6.2, the action on the line L^{123}_{mn} is the same as
\\ the action on the pair (t_1, t_3/t_2).
\\ The matrices returned by the following few functions are permutation
\\ matrices.
\\ Represent the element g on the basis { e^n t_j/t_i | n=1,3,5,7 }
\\ Here g is a 4-element vector representing an element of Gamma', and
\\ the return value is given by Lemma 3.20. And t_0 means 1.
\\ For our example, we have:
\\ g_1^4 fixes everything.
\\ g_2^2 fixes e^n t_1, and takes e^n t_3/t_2 -> e^(n+4) t_3/t_2
\\ takes e^n t_2 -> e^(n+4) t_2, fixes e^n t_1/t_3
\\ fixes e^n t_3, takes e^n t_2/t_1 -> e^(n+4) t_2/t_1
\\ g_2.g_3 fixes e^n t_1 and e^n t_3/t_2
\\ takes e^n t_2 -> e^(n+2) t_2, e^n t_1/t_3 -> e^(n-2) t_1/t_3
\\ takes e^n t_3 -> e^(n+2) t_3, e^n t_2/t_1 -> e^(n+2) t_2/t_1
\\ g_2.g_4 takes e^n t_1 -> e^(n+4) t_1, e^n t_3/t_2 -> e^(n+2) t_3/t_2
\\ takes e^n t_2 -> e^(n-2) t_2, e^n t_1/t_3 -> e^(n+4) t_1/t_3
\\ takes e^n t_3 -> e^(n+4) t_3, e^n t_2/t_1 -> e^(n-2) t_2/t_1
small_rep(g,i,j) =
{
local(b);
b = 2*g[4];
if(j > 0, b += g[j]);
if(i > 0, b -= g[i]);
b %= 4;
if(b == 0, matid(4),
if(b == 1, [0,0,0,1;1,0,0,0;0,1,0,0;0,0,1,0],
if(b == 2, [0,0,1,0;0,0,0,1;1,0,0,0;0,1,0,0],
[0,1,0,0;0,0,1,0;0,0,0,1;1,0,0,0]
)
)
);
}
\\ This function does the same for the element g_5.g, where g_5 is the
\\ fifth generator of Gamma which takes i -> -i.
\\
\\ Example: g_5.g_2.g_3 takes
\\ e^n t_1 -> e^(-n) t_1, e^n t_3/t_2 -> e^(-n) t_3/t_2
\\ e^n t_2 -> e^-(n+2) t_2, e^n t_1/t_3 -> e^-(n-2) t_1/t_3
\\ e^n t_3 -> e^-(n+2) t_3, e^n t_2/t_1 -> e^-(n+2) t_2/t_1
small_rep_conj(g,i,j) = [0,0,0,1;0,0,1,0;0,1,0,0;1,0,0,0] * small_rep(g,i,j);
\\ This function implements Algorithm 3.21.
\\ Represent the element g on the basis { L^{pqr}_{mn} | m,n=1,3,5,7 }
\\ Example: using the descriptions above, we get
\\ g_1^4 fixes everything
\\ g_2^2 : L^{123}_{mn} -> L^{123}_{m,n+4}
\\ L^{231}_{mn} -> L^{231}_{m+4,n}
\\ L^{312}_{mn} -> L^{312}_{m,n+4}
\\ g_2.g_3 : fixes L^{123}_{mn}
\\ L^{231}_{mn} -> L^{231}_{m+2,n-2}
\\ L^{312}_{mn} -> L^{312}_{m+2,n+2}
\\ g_2.g_4 : L^{123}_{mn} -> L^{123}_{m+4,n+2}
\\ L^{231}_{mn} -> L^{231}_{m-2,n+4}
\\ L^{312}_{mn} -> L^{312}_{m+4,n-2}
big_rep(g,p) = tensor(small_rep(g,(p%3)+1,((p+1)%3)+1), small_rep(g,0,p));
\\ ... and the same for g_5.g
\\ g_5.g_2.g_3 : L^{123}_{mn} -> L^{123}_{-m,-n}
\\ L^{231}_{mn} -> L^{231}_{-(m+2),-(n-2)}
\\ L^{312}_{mn} -> L^{312}_{-(m+2),-(n+2)}
big_rep_conj(g,p) = tensor(small_rep_conj(g,(p%3)+1,((p+1)%3)+1), \
small_rep_conj(g,0,p));
\\ Put three 16x16 matrices together into a 48x48 one
concat48(M1,M2,M3) =
{
M1 = concat(M1~,matrix(matsize(M1)[2],32))~;
M2 = concat(matrix(matsize(M2)[2],16),
concat(M2~,matrix(matsize(M2)[2],16)))~;
M3 = concat(matrix(matsize(M3)[2],32),M3~)~;
concat(M1,concat(M2,M3));
}
\\ Now we can get the representation of the element g on the whole
\\ group generated by the 48 lines. The result of this is a great big
\\ 48x48 permutation matrix.
rep48(g) = concat48(big_rep(g,1), big_rep(g,2), big_rep(g,3));
rep48_conj(g) = concat48(big_rep_conj(g,1), big_rep_conj(g,2), \
big_rep_conj(g,3));
\\ Represent a group element on NS(V)
\\ Here we use the basis chosen at the beginning of this file for NS(V).
\\ The matrix E is the quotient map onto NS(V); the matrix D is a section
\\ of it, that is, a right inverse.
ns_rep(g) = E*rep48(g)*D;
ns_rep_conj(g) = E*rep48_conj(g)*D;
\\ Represent a group on NS(V).
\\ This function puts together the representations of each generator
\\ of the subgroup. It returns a vector of 20x20 matrices, one matrix
\\ for each generator. This is what's needed for the algorithms found
\\ in h1.gp.
whole_ns_rep(H) =
{
local(G, n);
n = matsize(H[1])[2];
G = vector(n, i, ns_rep(H[1][,i]));
if(H[2] != [], G = concat(G, [ns_rep_conj(H[2])]));
G;
}
\\ Does H fix F?
fixes(H,F) =
{
local(G,n);
G = whole_ns_rep(H);
n = matsize(G)[2];
for(i=1,n,if(G[i]*F != F, return(0)));
1;
}
\\ Find a basis for H^0(H, NS(V))
\\ In our example, we find that NS(V) is generated by the classes of
\\ a plane section and of
\\ F = L^{123}_{11} + L^{123}_{15} + L^{123}_{31} + L^{123}_{35}
\\ which is the class of a fibration in curves of genus 1.
ns(H) = H0(whole_ns_rep(H));
\\ Find a set of generators for d^-1(H^1(H, NS(V)))
\\ See h1.gp for more explanation
\\ For our example, this shows that H^1(H, NS(V)) is of order 2, and
\\ gives the generator.
H1(H) =
{
local(X);
X = H1torsion(whole_ns_rep(H), 256); X[2] /= 256;
X;
}
\\ ---------------------------------------------------------------------
\\ Next we turn to finding the subgroup H corresponding to a surface,
\\ defined by coefficients a_i. This is described in Section 3.4 of
\\ the thesis.
\\ Take a factorisation and turn it into a vector on the basis S,
\\ where S is a vector of primes in ascending order.
factor2basis(f, S) =
{
local(j, n, v);
j = 1;
n = matsize(f)[1];
v = vector(matsize(S)[2]);
for(i = 1, n,
while(S[j] < f[i,1], j++);
v[j] = f[i,2];
);
v~;
}
\\ This is Algorithm 3.10. The reasoning behing the algorithm is not
\\ trivial and can be found in Section 3.4. Vectors and matrices will
\\ represent elements of Q*/Q*^4, or more specifically of that part of
\\ it generated by the primes in S.
\\ Our running example has [a0, a1, a2, a3] = [1, 1, -6, -12].
\\ In the comments for this function, 4[x] indicated the fourth root of x.
find_group(a0, a1, a2, a3) =
{
local(f1, f2, f3, S, n, Q, Qe, K1, K2, K3, v1, v2, v3, M, v, w1, w2, w3);
f1 = factor(a1/a0);
f2 = factor(a2/a0);
f3 = factor(a3/a0);
\\ First find out which primes interest us
S = Set([-1, 2]);
S = setunion(S, Set(f1[,1]));
S = setunion(S, Set(f2[,1]));
S = setunion(S, Set(f3[,1]));
n = matsize(S)[2];
for(i = 1, n, S[i] = eval(S[i]));
\\ Sort them - you might think that has been done already, but
\\ at the moment they're in lexicographic order rather than numerical
\\ order. This leads to a rather obscure bug.
S = vecsort(S);
\\ In our example, we have S = [-1, 2, 3].
\\ Now find the subgroups of Q* containing fourth powers of
\\ the various fields
Q = 4 * matid(n); Q[1,1] = 2; \\ In Q, all fourth powers are obvious
\\ [2 0 0]
\\ Q = [0 4 0] representing 1, 16, 81
\\ [0 0 4]
\\ In Q(e), we find that 4 and -4 are both fourth powers.
Qe = Q; Qe[1,2] = 1; Qe[2,2] = 2; Qe[2,1] = 2;
\\ [2 1 0]
\\ Qe = [2 2 0] representing 4, -4, 81
\\ [0 0 4]
\\ Now Qe is a basis for (Q(e)*^4 intersect Q*) / Q*^4
\\ After this, Kummer theory tells us that the new fourth powers
\\ in each extension are just the obvious ones.
v3 = factor2basis(f3, S); K3 = concat(v3, Qe);
\\ Now K3 is a set of generators for
\\ (Q(e, 4[-12])*^4 intersect Q*) / Q*^4
\\ So in fact we have
\\ [1 2 1 0]
\\ K3 = [2 2 2 0] representing -12, 4, -4, 81.
\\ [1 0 0 4]
v2 = factor2basis(f2, S); K2 = concat(v2, K3);
\\ [1 1 2 1 0]
\\ K2 = [1 2 2 2 0] representing -6, -12, 4, -4, 81.
\\ [1 1 0 0 4]
v1 = factor2basis(f1, S); K1 = concat(v1, K2);
\\ [0 1 1 2 1 0]
\\ K1 = [0 1 2 2 2 0] representing 1, -6, -12, 4, -4, 81.
\\ [0 1 1 0 0 4]
\\ Of course these matrices have many redundant columns.
\\ K1 and K2 are both equivalent to the identity, meaning that the
\\ corresponding extensions contain fourth roots of both 2 and 3.
\\ K3 is equivalent to diag[1,2,1], indicating that the extension
\\ contains a fourth root of 3 and a square root of 2.
\\ Work out the degrees of the relevant extensions
M = matrix(4,4); v = vector(4);
v[1] = matdet(mathnf(K1));
v[2] = matdet(mathnf(K2));
v[3] = matdet(mathnf(K3));
v[4] = matdet(Qe);
M[1,1] = 4 * v[1] / v[2];
M[2,2] = 4 * v[2] / v[3];
M[3,3] = 4 * v[3] / v[4];
M[4,4] = 1;
\\ Our example has [4 0 0 0] at this stage.
\\ M = [0 2 0 0]
\\ [0 0 1 0]
\\ [0 0 0 1]
\\ At this point we don't know what the elements above the diagonal are.
\\ We write H = Gal(Q(e,4[-12],4[-6],4[1]) / Q).
\\ This matrix means that:
\\ Gal(Q(e,4[-12],4[-6],4[1]) / Q(e,4[-12],4[-6])) is generated
\\ by an automorphism taking 4[1] -> 4[1] (that is, the identity);
\\ Gal(Q(e,4[-12],4[-6]) / Q(e,4[-12])) is generated by an
\\ automorphism taking 4[-6] -> -(4[-6]) and whose lift to H has
\\ some action, as yet unknown, on 4[1];
\\ Gal(Q(e,4[-12]) / Q(e)) is generated by an automorphism taking
\\ 4[-12] -> i(4[-12]) and whose lift to H has some unknown action
\\ on 4[-6] and 4[1];
\\ Gal(Q(e)/Q(i)) is generated by an automorphism taking \sqrt{2}
\\ to -\sqrt{2} and whose lift to H has some action on all the other
\\ fourth roots around.
\\ We now need to fill in the above-diagonal elements.
\\ Find out how the various extensions are related to each other
w1 = matsolvemod(K2, 4, (4 / M[1,1]) * v1);
\\ In our example, this tells us that 4[1] is 1. (w1 = [0,0,0,0,0]~)
w2 = matsolvemod(K3, 4, (4 / M[2,2]) * v2);
\\ We already know from before that 4[-6]^2 lies in Q(e,4[-12]). Here
\\ we find that 4[-6]^2 = 4[-12]^2 / (4[4] * 4[-1]^2) (*)
\\ so w2 = [2,-1,-2,0]~
w3 = matsolvemod(Qe, 4, (4 / M[3,3]) * v3);
\\ And here we find that -12 lies in Q(e) (no surprise there).
\\ So w3 = [0,0,0]~.
\\ Work out the Galois group
M[1,2] = M[2,2] * w1[1] * M[1,1] / 4;
\\ So the automorphism taking 4[-6] -> -4[-6] has to fix 4[1]. Good.
M[2,3] = M[3,3] * w2[1] * M[2,2] / 4;
\\ So the automorphism taking 4[-12] -> i(4[-12]) must take 4[-6] -> -4[-6],
\\ by looking at the expression (*) above
M[1,3] = (M[2,3] * w1[1] + M[3,3] * w1[2]) * M[1,1] / 4;
\\ and also fix 4[1].
M[3,4] = 2 * w3[1] * M[3,3] / 4;
\\ Finally, the automorphism taking \sqrt{2} -> -\sqrt{2} may fix 4[-12]
M[2,4] = (M[3,4] * w2[1] + 2 * w2[2]) * M[2,2] / 4;
\\ and take 4[-6] -> -i(4[-6]) (from (*))
M[1,4] = (M[2,4] * w1[1] + M[3,4] * w1[2] + 2 * w1[3]) * M[1,1] / 4;
\\ and fix 4[1].
\\ Lastly, we have a fifth generator which takes i -> -i. As we're only
\\ interested in this element up to conjugacy, that means that we only
\\ need to know its entries modulo 2. It just depends on the sign
\\ of the a_i/a_0, which we have already found earlier.
v[1] = v1[1] % 2;
v[2] = v2[1] % 2;
v[3] = v3[1] % 2;
v[4] = 0;
\\ Our example has v = [0,1,1,0].
[mathnf(M), v];
}
\\ ---------------------------------------------------------------------
\\ M = intersection matrix for NS(V), so non-singular and of determinant
\\ minus a power of 2.
\\ This works out whether V contains a pencil of curves of genus 1, that
\\ is, whether M represents 0 non-trivially. See pp43--44 of Serre for more
\\ info and proofs.
hasgenus1(M) =
{
\\ rank >= 3: yes
if(matsize(M)[1] >= 3, 1,
\\ rank <= 1: no
if(matsize(M)[1] == 1, 0,
\\ rank = 2: if -d is a square
if(matsize(M)[1] == 2, issquare(-matdet(M)))
)
)
}
\\ Return the class of a genus 1 curve on V. This only works when
\\ dim NS(V) = 2, as higher dimensions involve solving higher-dimensional
\\ quadratic forms, which we don't want to get into.
\\ N = basis for NS(V), in terms of D.
\\ For our example, this finds the fibre F given above. Or maybe 2P-F,
\\ where P is a plane section; that works as well.
findgenus1(N) =
{
local(M, d, n);
M = N~*D~*Q*D*N;
M /= content(M);
if(matsize(N)[2] <= 1, [],
if(matsize(N)[2] == 2,
d = matdet(M);
if(issquare(-d,&n), N*[-M[2,1]+n, M[1,1]]~, [])
,
[]
);
);
}
\\ This is Algorithm 2.20.
\\ Return all divisor classes of degree d and self-intersection at most n
listdivisors(N,d,n) =
{
local(T,U,M,v,L,s,r,t);
if(4*n > d*d, [],
r = matsize(N)[2];
M = N~*D~*Q*D*N;
for(i = 1, r, if(M[1,i] != 0, s = i));
if(s > 1,
T = matrix(r, r-1, i, j, if(i == 1, -M[1,j+1]/4, if(i == j+1, 1, 0)));
U = matid(r-1); U[s-1,] = 16 * T[1,] / M[1,s]; U[s-1,s-1] = 4 / M[1,s];
v = vectorv(r-1); v[s-1] = d/4;
M = U~*T~*M*T*U;
t = fincke_pohsti(-M, d*d/4, v, 1)[1];
L = T*U*fincke_pohsti(-M, d*d/4, v, t)[3];
);
L
);
}