"""
The purpose of this program is to compute the Weyl group
of a cubic surface defined over Q using the algorithm
of Andreas-Stephan Eisenhans and Joerg Jahnel:
"Experiments with general cubic surfaces."
Usage: sage wg.sage -r 100 - b 20
This command applies the Elsenhans-Jahnel algorithm
to 20 cubic surfaces
x^3 + y^3 + z^3 + w^3 + ayzw + bxzw + cxyw + dxyz
with a, b, c, d in the range -100...100. In the
reported results, "W" means that the Galois
group is the full Weyl group of E6.
Also:
sage wg.sage -c 1 2 3 4
- test the cubic surface with parameters 1 2 3 4
sage wg.sage -h
- display help
"""
from random import random
import sys
from optparse import OptionParser
usage = "%prog options FILES"
parser=OptionParser(usage)
parser.add_option("--random", "-r",
dest="random",
help="sage wg.sage -r 10 -c 5: test 10 random cubic surfaces with |coeffs| <= 5")
parser.add_option("--coeffbound", "-b",
dest="coeffbound",
help="sage wg.sage -r 10 -c 5: test 10 random cubic surfaces with |coeffs| <= 5")
parser.add_option("--coefficients", "-c",
dest="coefficients",
action = "store_true",
help="sage wg.sage -c 2 3 5 7 tests cubic surface with parameters 2, 3, 5, 7")
parser.add_option("--verbose", "-v",
dest="verbose",
action="store_true",
help="verbose output")
parser.add_option("--info", "-i",
dest="info",
action="store_true",
help="print info on meaning of output")
parser.add_option("--nprimes", "-p",
dest="nprimes",
default="1000",
help="number of primes to use for testing (default = 1000)")
(options, args) = parser.parse_args()
###############################################
# Here are some utility functions for testing
# equality of lists and membership of a list in
# a set of lists.
###############################################
def listequal(A,B):
# return True iff A == B as sets
if len(A) != len(B):
return False
else:
A.sort()
B.sort()
same = True
for k in range(0,len(A)):
if A[k] != B[k]:
same = False
return same
def listinset(L,S):
# L is a list and S is a set of lists.
# Return True iff L is an element of S
isinset = False
for k in range(0,len(S)):
if listequal(L,S[k]):
isinset = True
return isinset
A = [[9,9,9]]
B = [[1,1,5,5,5,5,5],[2,5,5,5,10]]
C = [[1,4,4,6,12],[2,5,5,5,10],[1,2,8,8,8]]
types = [ [A,"A"], [B,"B"], [C,"C"]]
def type(L):
# return string "A", "B", "C"
# according to whether the list of degrees
# L is in the set of lists A, B, or C
theTypes = ""
for T in types:
if listinset(L,T[0]):
theTypes = theTypes + T[1]
return theTypes
def update_degrees(degrees, typestring):
# Add type of current list of degrees to typestring
theType = type(degrees)
if "A" not in typestring and "A" in theType:
typestring += "A"
if "B" not in typestring and "B" in theType:
typestring += "B"
if "C" not in typestring and "C" in theType:
typestring += "C"
return typestring
def multiplicities(F):
# return list of multiplicities of a factorization F
return map( lambda x: x[1], F)
######################################
# main function follows. It returns
# a string which gives both intermediate
# and final information on the Elsenhans-Jahnel
# algorithm. If the string contains "W", then
# the image of the Galois group is the Weyl
# group of E6.
def wg(aa,bb,cc,dd):
# Set up two polynomial rings, one in variables t, a, b, c, d
# the other in variables x, y, z, w
# orderings are 'lex', 'revlex', 'degrevlex', etc.
R. = PolynomialRing(QQ, 5, order = 'lex')
S. = PolynomialRing(QQ, 4)
# The cubic polynomial we will use:
f = x^3 + y^3 + z^3 + w^3 + aa*y*z*w + bb*x*z*w + cc*x*y*w + dd*x*y*z
# print f # debug
# l(t) is a general line.
# The parameters a,b,c,d are parameters in
# the Grassmannian G(2,4)
l = lambda t: (1,t,a+b*t,c+d*t)
# fl(t) is the composition of f(x,y,z,w) and l(t)
fl = lambda t: f(l(t))
# I is the ideal generated by fl(0), fl(1), fl(2), fl(3)
I = ideal(fl(0),fl(1),fl(2),fl(3))
# B is the Groebner basis of I
B = I.groebner_basis();
g = B[0]
"""
Now comes the main code. We start with an initial prime p. Then we
reduce the polynomial g modulo p, coercing it into a polynomial ring
in one variable at the same time. NOTE: we assume the outcome of the
computation gives a univariate polynomial in d. Then we find the factors of h,
which we store in fh. From fh we compute a list of degrees of the factors.
According to the Eisenhans-Jahnel paper, the list of degrees is of type A, B, C
of of undefined type. The type is computed for the given prime. Finally,
we increment the prime as long as required.
"""
typestring = ""
if g.degree() < 27:
return "S"
else: # degree is 27, so continue
p = 101
N = 0 # N is the index current test, e.g., we reduce modulo the Nth prime
NPRIMES = int(options.nprimes) # number of primes to use
while ("W" not in typestring) and (N < int(options.nprimes)):
# "W" signals success: we know that the image of the
# Galois group is E6
N = N + 1
p = next_prime(p)
L = GF(p)
T.__ = PolynomialRing(QQ)
# Coerce g to ZZ[u] and reduce g modulo p:
phi = R.hom([0,0,0,0,u], T)
h = phi(g) # h is in QQ[u]
hc = h.coeffs()
hd = map( lambda x: x.denominator(), hc)
hlcmd = lcm(hd)
h = h*hlcmd # now h is in ZZ[u]
if h.leading_coefficient() % p == 0:
typestring += "P"
else: # p does not divide leading coefficient so continue
fh = factor(h.change_ring(L))
if max(multiplicities(fh)) > 1:
typestring += "M"
else: # no multiple factors, so continue
degrees = map( lambda x: x[0].degree(), fh)
typestring = update_degrees(degrees, typestring)
# print p, degrees, type(degrees), typestring # debug
if ("A" in typestring) and ("B" in typestring) and ("C" in typestring):
typestring += "W"
typestring = `N`+":"+typestring
return typestring
################################################################################
def info():
print
print " Is image of Gal(bar Q/Q) in NS(Cubic surface) = W(E_6)?"
print " Equation is x^3 + y^3 + z^3 + w^3 + ayzw + bxzw + cxyw + dxyz"
print " where a, b, c, d are displayed below"
print
print " Codes:"
print " S - degree of polynomial < 27 (fatal)"
print " P - leading coefficient was divisible by p (not fatal)"
print " M - there were multiple factors mod p (not fatal)"
print " ABC - one of these types occured (good)"
print " 'W' means that imaage of Galois group is W(E_6)"
print " '' means that the computations don't prove anything"
print
def run( MAX, N ):
rand = lambda : int(2*random()*MAX - MAX)
hits = 0
if options.verbose:
info()
for i in range(0,N):
a,b,c,d = rand(), rand(), rand(), rand()
result = wg(a,b,c,d)
if "W" in result:
status = "E6"
hits += 1
else:
status = ""
if options.verbose:
print "%5d. %5d %5d %5d %5d %15s %2s" % (i+1, a,b,c,d, result, status)
else:
print "%5d. %5d %5d %5d %5d: %2s" % (i+1, a, b, c, d, status)
print; print "%5d hits out of %d" % (hits, N)
if options.random and options.coeffbound:
run( int(options.coeffbound), int(options.random) )
if options.coefficients:
a, b, c, d = map( lambda x: int(x), args )
result = wg(a,b,c,d)
if options.verbose:
print " ", result
else:
if "W" in result:
print " G = E6"
else:
print " ???"
if options.info:
info()__