"""
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.<t,a,b,c,d> = PolynomialRing(QQ, 5, order = 'lex')
  S.<x,y,z,w> = 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.<u> = 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()