# Constants
#
#  Q8: The Quaternion group.
#  a,b,a2,a3,ab,a2b,a3b: Quaternion elements (short form of a^2*b, ...). 
#
# Functions:
#
#  c2G( v ): makes new G element
#  G2str(x): prints G elements
#  Gu(): Element of order 2
#  commutator(x,y): commutator
#  swapper(x,y): swapper
#  wt(bl,x): weight
#  T(C): Subgroup of order 2
#  K(C): Kernel
#  S(C): Lineal spawn
#  c_type(C): type of C
#  is_hadamard(C,bl): verifies if Hadamard code
#  shape(C): shape of C
#  analysis(C,bl): utility, main C characteristics
#
# Basic example:
#  load "newlib.sage"
#
#  G=direct_product_permgroups( [Q8]*8 )
#  bl=[8]*8
#
#  r1=c2G( [ 1,a,a2,a3,1,a,a2,a3 ] );
#  s1=c2G( [ 1,b,a2,b,a2,ab,1,ab ] );
#  C=G.subgroup([r1,s1])
#  analysis(C,bl)
  


# Quaternonic
Q8=DiCyclicGroup(2)
cn=8

# DihedralGroup
#Q8=DihedralGroup(4)
#cn=4

a=Q8.gen(0)
b=Q8.gen(1)
a2=a^2
a3=a^3
ab=a*b
a2b=a^2*b
a3b=a^3*b

#
# core functions
#

def tuple_in(mi,x):
  for v in list(x):
    if v<mi or v>=mi+cn: return False
  return True 

def tuple_down(mi,x):
  return tuple([ i-mi for i in list(x) ]) 

def tuple_up(mi,x):
  return [ i+mi for i in list(x) ]

def cycle2component( v, mi ):
  ct = v.cycle_tuples()
  return Q8([ tuple_down(mi,x) for x in ct if tuple_in(mi+1,x) ])

#
# I/O functions
#

def c2G( v ):
  r"""
  Converts a list of components to a G group element

  INPUT:
  - v -- Lists of components

  OUTPUT:

  G group element

  EXAMPLES:

  G=direct_product_permgroups( [Q8]*8 )
  x=c2G( [ 1, a, a^2, a^3, b, a*b, a^2*b, a^3*b ] )
  G2str(x)
  y=c2G( [ 1, a, a2, a3, b, ab, a2b, a3b ] )
  G2str(y)
  """

  r=[]
  mi=0
  for i in v:
    if i == 1:  
      i=a^0
    ct=i.cycle_tuples()
    r.extend( [ tuple([ n+mi for n in list(t) ]) for t in ct ] )
    mi += cn
  return G(r)

def G2str_aux2(x):
#  if isinstance(x,list): return [ G2str(i) for i in x]
  if x==a^0: return "1"
  if x==a^1: return "a"
  if x==a^2: return "a^2"
  if x==a^3: return "a^3"
  if x==b: return "b"
  if x==a*b: return "a*b"
  if x==a^2*b: return "a^2*b"
  if x==a^3*b: return "a^3*b"
  return "??"

def G2str_aux(v):
  if v==None: return 'None'
  if isinstance(v,list): 
    return [ G2str(i) for i in v]
  if isinstance(v,set): 
    return [ G2str(i) for i in v]
  n=v.parent().degree()/cn
  r=range(0,n*cn,cn)
  res=[ G2str_aux2(cycle2component(v,mi)) for mi in r ]
  return res

def G2str(x):
  r"""
  Return representation of a G element (Latex like).

  INPUT:
  - x ``x`` -- Group element(s).

  OUTPUT:

  Human readable representation.

  EXAMPLES:

  G=direct_product_permgroups( [Q8]*8 )
  x=c2G( [ 1, a, a^2, a^3, b, a*b, a^2*b, a^3*b ] )
  G2str(x)
  """

  res=G2str_aux(x)
  return "[ "+(", ".join(res))+" ]"
  
#
# basic concepts 
#

def Gu():
  r"""
  Return element of order two of G.

  INPUT:

  - x -- Group element.
  - y -- Group element.

  OUTPUT:

  The element of order 2 in G

  EXAMPLES:

  G=direct_product_permgroups( [Q8]*8 )
  u=Gu()
  G2str(u)
  """

  n=G.degree()/cn
  res=[]
  ct=(a^2).cycle_tuples()
  for mi in range(0,n*cn,cn):
    res.extend( [ tuple([ x+mi for x in list(t) ]) for t in ct ] )
  return G(res)

def commutator(x,y):
  r"""
  Return the commutator of two elements, [x,y].

  INPUT:

  - x ``x`` -- Group element.
  - y ``y`` -- Group element.

  OUTPUT:

  The commutator [x,y].

  EXAMPLES:

  G=direct_product_permgroups( [Q8]*8 )
  x=c2G( [ 1, a, a^2, a^3, b, a*b, a^2*b, a^3*b ] )
  y=c2G( [ b, a*b, a^2*b, a^3*b, b, a*b, a^2*b, a^3*b ] )
  G2str( commutator(x,y) )
  """
  return x^-1*y^-1*x*y

def swapper_aux(x,y):
  if x==a^0 or x==a^2: return a^0
  if y==a^0 or y==a^2: return a^0
  if (x==a or x==a^3) and (y==a*b or y==a^3*b): return a^0
  if (x==b or x==a^2*b) and (y==a or y==a^3): return a^0
  if (x==a*b or x==a^3*b) and (y==b or y==a^2*b): return a^0
  return a*a

def swapper(x,y):
  r"""
  Return the commutator of two elements, [x,y].

  INPUT:

  - x ``x`` -- Group element.
  - y ``y`` -- Group element.

  OUTPUT:

  The commutator [x,y].

  EXAMPLES:

  G=direct_product_permgroups( [Q8]*8 )
  x=c2G( [ 1, a, a^2, a^3, b, a*b, a^2*b, a^3*b ] )
  y=c2G( [ b, a*b, a^2*b, a^3*b, b, a*b, a^2*b, a^3*b ] )
  G2str( swapper(x,y) )
  """

  n=x.parent().degree()/cn
  r=range(0,n*cn,cn)
  res = [ swapper_aux( cycle2component(x,mi), 
                       cycle2component(y,mi) ) for mi in r ] 
  return c2G( res )

def wt_aux(b,x):
  if b==2 and x==a^0: return 0
  if b==2: return 1
  if b==4 and x==a^0: return 0
  if b==4 and x==a^2: return 2
  if b==4: return 1
  if b==8 and x==a^0: return 0
  if b==8 and x==a^2: return 4
  if b==8: return 2
  return -1

def wt(bl,x):
  r"""
  Return the weight in bits of the group element x, 
  using component map bl.

  INPUT:

  - bl ``bl`` -- Component map.
  - x ``x`` -- Group element.

  OUTPUT:

  Integer, the weight of the element.

  EXAMPLES:

  G=direct_product_permgroups( [Q8]*8 )
  bl=[8]*16
  x=c2G( [ 1, a, a^2, a^3, b, a*b, a^2*b, a^3*b ] )
  wt(bl,x)
  """

  n=x.parent().degree()/cn
  res=0
  for i in range(0,n):
    mi = i*cn
    xi = cycle2component(x,mi)
    res += wt_aux(bl[i],xi)
  return res

#
# main subgroups  
#

def T(C):
  r"""
  The subgroup of all order 2 elements of C.

  INPUT:

  - C ``C`` -- Subgroup of G.

  OUTPUT:

  Subgroup of C with all order 2 elements.

  EXAMPLES:

  G=direct_product_permgroups( [Q8]*16 )
  r1=c2G( [ a ] *16 )
  r2=c2G( [ 1,a,a^2,a^3 ] *4 )
  C=G.subgroup([r1,r2])
  TC=T(C)
  len(TC) # =4
  """

  elems = [ x for x in C.list() if x^2==C.one() ]
  TMP=C.subgroup(elems)
  return C.subgroup(TMP.gens_small())

def K_aux(C,T,gens,x):
  for i in gens:
    if not(swapper(i,x) in T): return False
#    if not(swapper(i,x) in T): return False
#    if not(swapper(x,i) in T): return False
  return True

def K(C):
  r"""
  The kernel of C.

  INPUT:

  - C ``C`` -- Subgroup of G.

  OUTPUT:

  Kernel of C

  EXAMPLES:

  G=direct_product_permgroups( [Q8]*16 )
  r1=c2G( [ a ] *16 )
  r2=c2G( [ 1,a,a^2,a^3 ] *4 )
  C=G.subgroup([r1,r2])
  KC=K(C)
  len(KC) # = 16
  """

  TC=T(C)
  gens=C.gens_small()
  elems = [ x for x in C if K_aux(C,TC,gens,x) ]
  TMP=C.subgroup(elems)
  return C.subgroup(TMP.gens_small())

def S(C):
  r"""
  The lineal spawn of C.

  INPUT:

  - C ``C`` -- Subgroup of G.

  OUTPUT:

  Lineal spawn of C

  EXAMPLES:

  G=direct_product_permgroups( [Q8]*16 )
  r1=c2G( [ a ] *16 )
  r2=c2G( [ 1,a,a^2,a^3 ] *4 )
  C=G.subgroup([r1,r2])
  SC=S(C)
  len(SC) # = 16
  """

  g=[]
  gens=C.gens_small()
  for a in gens:
    g.append(a)
    for b in gens:
      g.append(swapper(a,b))
  TMP=G.subgroup(g)
  return G.subgroup(TMP.gens_small())

def c_type(C):
  r"""
  The type of the subgroup C.

  INPUT:

  - C ``C`` -- Subgroup of G.

  OUTPUT:

  ((sigma,delta,rho),(T(C),Z(C)))

  EXAMPLES:

  G=direct_product_permgroups( [Q8]*16 )
  r1=c2G( [ a ] *16 )
  r2=c2G( [ 1,a,a^2,a^3 ] *4 )
  C=G.subgroup([r1,r2])
  SC=S(C)
  len(SC) # = 16
  """

  Z=C.center()
  TC=T(C)

  st=TC.order()
  sz=Z.order()
  sc=C.order()

  sigma=int(log(st,2))
  delta=int(log(sz,2))-sigma
  rho=int(log(sc,2))-delta-sigma

  return ((sigma,delta,rho),(TC,Z))
 
#
# Hadamard subgroups  
#

def is_hadamard(C,bl):
  r"""
  Verifies if C is a Hadamard code, 
  using component map bl.

  INPUT:

  - C ``C`` -- Subgroup of G.
  - bl ``Cl`` -- Component map

  OUTPUT:

  True if the code is Hadamard.
  If it is not a Hadamrd code, a non Hadamard element is returned

  EXAMPLES:

  G=direct_product_permgroups( [Q8]*16 )
  bl=[8]*16
  r1=c2G( [ a ] *16 )
  r2=c2G( [ 1,a,a^2,a^3 ] *4 )
  C=G.subgroup([r1,r2])
  is_hadamard(C,bl)
  """

  n=sum(bl)/2
  e0=False
  e1=False
  for x in C.list():
    wx=wt(bl,x)
    if wx==n:
      e1=True
      
    elif wx==0:
      e0=True

    elif wx != n/2: 
      return x
#  return e1 and e2
  return True

def mask_aux(x,y):
  if x==a^0 or x==a^2: return y
  else: return x

def mask(C,x,y):
  if commutator(x,y)!=C.one(): return None
  n=x.parent().degree()/cn
  r=range(0,n*cn,cn)
  res = [ mask_aux( cycle2component(x,mi), 
                    cycle2component(y,mi) ) for mi in r ] 
  return c2G( res )

def shape_sort(q1,m1,q2,m2):
  if q1==None and q2:
    return (q2,m2,None,m1)
  if q1==None and m1==None and m2:
    return (None,m2,None,None)
  return (q1,m1,q2,m2)

def shape_q(C,q,m,a):
  e=C.one()
  u=Gu()

  if q:                                          # X*
    co=commutator(q,a)
    if a^2==u: 
      if co==e: return (q,m,None)                   # => X*
      if co==u: return (q,m,a) # assert [q,a]==u    # => CONT a 
      return (q,m,q*a) # assert [q,a']==a'^2        # => CONT q*a

    if co!=e: return (q,m,a) # assert [q,a]==a^2    # => CONT a
    if m==None: return (q,a,None)                   # => XX
    return(q,mask(C,m,a),None)                      # => XX

  # q==None
  if m:                                          # _X
    nm=mask(C,m,a)
    if nm==None:
      return (None,m,a) # assert [m,a]!=e          # => CONT a
    if nm^2==u:
      return (nm,m,None)                           # => XX
    return (None,nm,None)                          # => _X

  # q==None and m==None                          # __
  if a^2==u: return (a,None,None)                  # => X_
  return (None,a,None)                             # => _X 

def shape_aux(C,q1,m1,q2,m2,a):
  e=C.one()
  u=Gu()

#  print("---IN---")
#  print(G2str(q1))
#  print(G2str(m1))
#  print(G2str(q2))
#  print(G2str(m2))
#  print(G2str(a))

  if a^2==e: return (q1,m1,q2,m2)

  (q1,m1,a) = shape_q(C,q1,m1,a)
  if a==None: return (q1,m1,q2,m2)
  (q2,m2,a) = shape_q(C,q2,m2,a)
  if a==None: return (q1,m1,q2,m2)
  (q1,m1,q2,m2) = shape_sort(q1,m1,q2,m2)

#  print("---POST---")
#  print(G2str(q1))
#  print(G2str(m1))
#  print(G2str(q2))
#  print(G2str(m2))
#  print(G2str(a))
#  print

  if a^2==u:
    if q2: return (q1,m1,q2,m2)              # X? X? => IGN
    if m2: return shape_aux(C,q1,m1,a,None,m2) # ?X _X => REDO ?X X_ % m2
    return (q1,m1,a,None)                    # ?? __ => ?? X_ 

  if m2==None:
    if q2==None: return (q1,m1,None,a)      # ?? __ => ?? _X 
    return (q1,m1,q1*q2,a)                  # X? X_ => X? XX

  return (q1,m1,q2,m2)                      # => IGN

def shape(C):
  r"""
  The shape of a code Hadamard code C.

  INPUT:

  - C ``C`` -- Subgroup of G.

  OUTPUT:

  (shape,eps) where:
  - shape: the shape of the C code.
  - eps: epsilon value.

  EXAMPLES:

  G=direct_product_permgroups( [Q8]*16 )
  r1=c2G( [ a ] *16 )
  r2=c2G( [ 1,a,a^2,a^3 ] *4 )
  C=G.subgroup([r1,r2])
  shape(C)
  """

  u=Gu()
  (q1,m1,q2,m2)=(None,None,None,None)

  for a in C.gens():
    st=shape_aux(C,q1,m1,q2,m2,a)
    if st==None: return None
    (q1,m1,q2,m2)=st

  (q1,m1,q2,m2) = shape_sort(q1,m1,q2,m2)

#  print(G2str(q1))
#  print(G2str(m1))
#  print(G2str(q2))
#  print(G2str(m2))

  if q1==None: 
    if m1==None:                 # (_,_) (_,_)
      return (1,0)
    if m2==None:                 # (_,X) (_,_)
      return (1,0) 
                                 # (_,X) (_,X)
    if m1^2==m2^2:
      return (4,1)
    return (-1,-1)                 

  if q2==None: # 
    if m1==None and m2==None:    # (X,_) (_,_)
      return (1,0)
    if m1==None and m2!=None:    # (X,_) (_,X)
      return (3,0)
    if m2==None:                 # (X,X) (_,_)
      return (1,0)
                                 # (X,X) (_,X)
    if m1^2==m2^2 or (m1*q1)^2==m2^2:
      return (4,1)
    return (-1,-1) 

  if m1==None:  
    return (2,1)                 # (X,_) (X,_)
  if m2==None:                 
    return (2,1)                 # (X,X) (X,_)
                                 # (X,X) (X,X)
  if m1^2==m2^2 or (m1*q1)^2==m2^2:
    return (5,2)  
  return (-1,-1) 

#
# tools 
#

def analysis(C,bl):
  r"""
  Prints the main characteristics of a subgroup C.

  INPUT:

  - C ``C`` -- Subgroup of G.
  - bl ``bl`` -- Component map.

  OUTPUT:

  None.

  EXAMPLES:

  G=direct_product_permgroups( [Q8]*16 )
  bl=[8]*16
  r1=c2G( [ a ] *16 )
  r2=c2G( [ 1,a,a^2,a^3 ] *4 )
  C=G.subgroup([r1,r2])
  analysis(C,bl)
  """

  lc=C.order()
  print "Dimension(C) = %d"%(log(lc,2))
  
  had=is_hadamard(C,bl)
  if had==True:
    print "Hadamard code."
  else:
    print "NOT a Hadamard code(example=%s)."%G2str(had)

  (ct,_)=c_type(C)
  print "type: sigma=%d,delta=%d,rho=%d"%ct
  sh=shape(C)
  print "shape=%d, eps=%d"%sh
  kc=K(C).order()
  print "kernel dimension = %d"%(log(kc,2))
  sc=S(C).order()
  print "rank = %d"%(log(sc,2))