###########################################
# Computation of the 40 Yoshida functions #
###########################################

# This file contains the computation of the 40 Yoshida functions. We use the date from the Root System E6 stored in "../../General" as well as the action of W(E6) on the ambient polynomial ring in the variables d1,...,d6.

# Since different versions of SAGE use a different ordering for the positive roots in the root system, we must load the sage objects for the two root systems stored in "../../General")

load("../../General/orbit.sage")

######################################################
# Setting up the ambient ring with the W(E6) action. #
######################################################

# The Yoshida functions will be elements of the following polynomial
# ring.
# A list of variables d_1, ..., d_6 for the coefficients
ds = [var("d%s"%i) for i in range(1,7)]
Sd =  PolynomialRing(QQ, ds)

#Sd.<d1,d2,d3,d4,d5,d6> = PolynomialRing(QQ, 6)

# Let D = <d1,...,d6>.

# Put the following bilinear form on D.
# (d_i, d_j) = 8/9 if i = j; -1/9 if i != j.

# Then the lattice spanned by the elements
# d_i - d_j (i,j distinct)
# +/- (d_i + d_j + d_k) (i,j,k distinct)
# +/- (d_1 + ... d_6)
# is isomorphic to the E6 lattice with the above 72 elements as the roots.

# We now identify the lattice above with the E6 lattice in sage.



######################################
# Data associated to the root system #
######################################

#We recall the relevant notation to have the simple roots alpha[i]'s, the positive roots, the Weyl group and its simple reflections for the root system E6.

RootSystemE6 = load("../../General/L6.sobj")

L6 = RootSystemE6.root_lattice()

# L6.dynkin_diagram()
#         O 2
#         |
#         |
# O---O---O---O---O
# 1   3   4   5   6
# E6
#diagram = [{1,3}, {3,4}, {4,2}, {4,5}, {5,6}]


# And the corresponding roots are labeled alpha[1], ..., alpha[6].

# We choose the isomorphism 

# alpha[1] -> d2 - d1
# alpha[2] -> d1 + d2 + d3
# alpha[3] -> d3 - d2
# alpha[4] -> d4 - d3
# alpha[5] -> d5 - d4
# alpha[6] -> d6 - d5

# This isomorphism is encoded in the following matrix

# Alpha to d matrix
AtoD = matrix([[-1, 1, 0, 0, 0, 0],
               [ 1, 1,-1, 0, 0, 0],
               [ 0, 1, 1,-1, 0, 0],
               [ 0, 0, 0, 1,-1, 0],
               [ 0, 0, 0, 0, 1,-1],
               [ 0, 0, 0, 0, 0, 1]])


# We write down the list of positive roots in the ds.

#pos = load("../../General/posRootsE6.sobj")
#positive_roots = [vector(ds).dot_product(AtoD*v.to_vector()) for v in pos]
# save(positive_roots, "Input/rootsE6.sobj")
positive_roots = load("Input/rootsE6.sobj")
roots = positive_roots

# For the record, here it is:
# [-d1 + d2,
#  d1 + d2 + d3,
#  -d2 + d3,
#  -d3 + d4,
#  -d4 + d5,
#  -d5 + d6,
#  -d2 + d4,
#  d1 + d2 + d4,
#  -d4 + d6,
#  -d1 + d3,
#  -d3 + d5,
#  -d2 + d5,
#  -d1 + d4,
#  d1 + d2 + d5,
#  -d3 + d6,
#  d1 + d3 + d4,
#  d2 + d3 + d4,
#  d1 + d2 + d6,
#  -d2 + d6,
#  -d1 + d5,
#  d1 + d3 + d5,
#  d1 + d4 + d5,
#  d2 + d3 + d5,
#  d1 + d3 + d6,
#  -d1 + d6,
#  d2 + d3 + d6,
#  d2 + d4 + d5,
#  d1 + d4 + d6,
#  d1 + d5 + d6,
#  d3 + d4 + d5,
#  d2 + d4 + d6,
#  d2 + d5 + d6,
#  d3 + d4 + d6,
#  d3 + d5 + d6,
#  d4 + d5 + d6,
#  d1 + d2 + d3 + d4 + d5 + d6]

##############
# Weyl Group #
##############

# The following are the generators of the Weyl group (in terms of
# alphas).

weyl_gens = L6.weyl_group().gens()

# We first convert them in terms of d's.
weyl_gens_ds = [AtoD * g * AtoD.inverse() for g in weyl_gens]

# In terms of the ds, all of them but one look like simple
# transpositions of the subscripts.

# The following helper function will let us convert these matrices
# into homomorphisms of Sd.

# Inputs:
#   matrix (a dxd matrix)
#   varlist (a list of d variables)

# Output:
#   a list of d linear combinations of the variables.
#   The i th element is the linear combination Sum matrix_{i,j} varlist_j

def varListFromMatrix(matrix, varlist):
    def newVar(col):
        return col * vector(varlist)
    return map(newVar, matrix.columns())

# These are the isomorphisms of Sd induced by the elements of weyl_gens_ds.
homs = {i:Sd.hom(varListFromMatrix(weyl_gens_ds[i], Sd.gens())) for i in range(0, len(weyl_gens_ds))}

#save(homs, "homOfDsFromWeylGroup.sobj")

# For the record, these are the results written explicitly

# {0: Ring endomorphism of Multivariate Polynomial Ring in d1, d2, d3, d4, d5, d6 over Rational Field
#   Defn: d1 |--> d2
#         d2 |--> d1
#         d3 |--> d3
#         d4 |--> d4
#         d5 |--> d5
#         d6 |--> d6,
#  1: Ring endomorphism of Multivariate Polynomial Ring in d1, d2, d3, d4, d5, d6 over Rational Field
#   Defn: d1 |--> 1/3*d1 - 2/3*d2 - 2/3*d3
#         d2 |--> -2/3*d1 + 1/3*d2 - 2/3*d3
#         d3 |--> -2/3*d1 - 2/3*d2 + 1/3*d3
#         d4 |--> 1/3*d1 + 1/3*d2 + 1/3*d3 + d4
#         d5 |--> 1/3*d1 + 1/3*d2 + 1/3*d3 + d5
#         d6 |--> 1/3*d1 + 1/3*d2 + 1/3*d3 + d6,
#  2: Ring endomorphism of Multivariate Polynomial Ring in d1, d2, d3, d4, d5, d6 over Rational Field
#   Defn: d1 |--> d1
#         d2 |--> d3
#         d3 |--> d2
#         d4 |--> d4
#         d5 |--> d5
#         d6 |--> d6,
#  3: Ring endomorphism of Multivariate Polynomial Ring in d1, d2, d3, d4, d5, d6 over Rational Field
#   Defn: d1 |--> d1
#         d2 |--> d2
#         d3 |--> d4
#         d4 |--> d3
#         d5 |--> d5
#         d6 |--> d6,
#  4: Ring endomorphism of Multivariate Polynomial Ring in d1, d2, d3, d4, d5, d6 over Rational Field
#   Defn: d1 |--> d1
#         d2 |--> d2
#         d3 |--> d3
#         d4 |--> d5
#         d5 |--> d4
#         d6 |--> d6,
#  5: Ring endomorphism of Multivariate Polynomial Ring in d1, d2, d3, d4, d5, d6 over Rational Field
#   Defn: d1 |--> d1
#         d2 |--> d2
#         d3 |--> d3
#         d4 |--> d4
#         d5 |--> d6
#         d6 |--> d5}




###################################
# Computing the Yoshida functions #
###################################

# We start with the one Yoshida function given by

# Yoshida1 = (d5 - d6) * (d4 + d5 + d6) * (-d3 + d4) * (d3 + d5 + d6) * (d2 + d3 + d4) * (-d1 + d2) * (d1 + d3 + d4) * (d1 + d2 + d6) * (d1 + d2 + d5)

# The rest are obtained as the W(E6) orbit, as follows:
# We use Yoshida1 and the homomorphisms from E6 to construct all 40 Yoshida functions (and their negatives). The file is stored as a SET and each Yoshida function is presented as a degree 9 polynomial in the polynomial ring S above.


#Yoshida_functions = computeOrbit(set([Sd(Yoshida1)]), homs)
# Takes about 1 minute on Quad core Intel(R) Core(TM) i5 with 4MB
# cache and 4GB RAM.

# The result is saved in "Yoshida_functions.sobj".
# This is a set of 80 functions (40 up to sign).  Each element is a degree 9 expanded polynomial in the variables d1, ..., d6 with coefficients in Q
#save(Yoshida_functions, "Yoshida_functions.sobj")


Yoshida_functions = load("../Output/Yoshida_functions.sobj")


################################
# Computing the Yoshida matrix #
################################

# The following function computes the 40x36 Yoshida_matrix 
# Convention: rows = exponents of the map: (roots of E6) --> Yoshida functions

# By definition, each Yoshida function is a product of roots.
# Up to sign, it is a product of positive roots.
# The Yoshida matrix encodes this product structure.

# We use the (ordered) list of positive roots stored in rootsE6.sobj

# The following function converts an exponent vector to a product of roots.
# We use it for getting Yoshida functions from the rows of the Yoshida matrix.
def VectorToYoshida(v, roots):
    return prod([roots[i]^v[i] for i in range(0,len(v))])


# The following helper function will let us compute the Yoshida matrix.

# Inputs:
#     f: A function
#     roots: A list of functions in the d's (the positive roots of E6)
# We assume that f factors as a product of elements of roots.

# Output:
#     A vector v such that f = Product_j (roots[j]^(v[j]))
def exponent_vector(f, roots):
    v = zero_vector(len(roots))
    for (r, i) in f.factor():
        if r in roots:
            root_location = roots.index(r)
        elif -r in roots:
            root_location = roots.index(-r)
        else:
            raise KeyError("Could not find a factor of f in roots: f = " + str(f) + ", roots = " + str(roots))

        v[root_location] =  v[root_location] + i

    return v

    
# We use the above to form the Yoshida matrix
# Note that repetitions (up to signs) are dropped.
# So the resulting matrix has 40 rows and 36 columns.

# exponent_vectors = []
# for f in Yoshida_functions:
#     v = exponent_vector(f, positive_roots)
#     if not v in exponent_vectors:
#         exponent_vectors.append(v)

# Yoshida_matrix = matrix(exponent_vectors)
# We save the result as a sage object
# save(Yoshida_matrix, "../Output/Yoshida_matrix.sobj")

# Caveat: Yoshida_functions does not have a natural ordering.
#         Therefore, Yoshida_matrix above may not agree on the nose
#         with the one stored in Yoshida_matrix.sobj. However, they
#         must have the same set of rows.

Yoshida_matrix = load("../Output/Yoshida_matrix.sobj")



# We check that the Yoshida_matrix has 9 ones per row and 10 ones per column, and that it contains only 0/1 entries:

# all([Yoshida_matrix[i][k] in [0,1] for i in range(0, Yoshida_matrix.nrows()) for k in range(0, Yoshida_matrix.ncols())])
# True
# all([sum(Yoshida_matrix[i]) == 9 for i in range(0, Yoshida_matrix.nrows())])
# True
# all([sum(Yoshida_matrix.transpose()[i]) == 10 for i in range(0, Yoshida_matrix.ncols())])
# True

# The Yoshida matrix has rank 16
# print Yoshida_matrix.rank()
# 16

# We check the index of the Z-span of the rows of the Yoshida_matrix  in Z^40. For this, it suffices to compute the Smith normal form of the Yoshida_matrix and check the entries of the diagonal.

# D, U, V = Yoshida_matrix.smith_form()
# print [D[i][i] for i in range(0,36) if D[i][i] != 0]
# [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3]
# This shows the index of the Z-span of the rows is 3.

# For the record, we print the Smith normal form of Yoshida_matrix:
# print D.str()
# [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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 3 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 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 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 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 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 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 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]
# [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 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 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 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 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 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 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 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 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]
# [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 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 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 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 0 0 0 0]




# In the future, we will always use the result stored in
# Yoshida_matrix.sobj as our Yoshida_matrix.


# By construction, Yoshidas[i] will always refer to the function defined by the i-th row of the Yoshida_matrix. In "OperatorsOnYoshidas.sage" this is encoded in the y_to_d variable.

# For the record, here is the matrix from "Yoshida_matrix.sobj":
# print Yoshida_matrix.str()
# [0 0 1 0 0 0 0 0 1 0 1 1 1 0 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 0 0 0 0 0 1]
# [0 1 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 1 0 0 0 0 0 1 1 1 0 0 0 0]
# [0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 1 0 1 1 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 0]
# [0 0 1 1 0 1 1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1]
# [0 1 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 0 0 0 0 1 0 1 0 1 0 0 0]
# [0 0 1 0 0 1 0 1 0 0 0 0 1 0 0 1 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 1 0]
# [1 0 0 1 0 1 0 0 0 0 0 0 0 1 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0]
# [0 0 0 0 0 0 1 0 1 1 1 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1]
# [0 0 0 0 0 1 1 1 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0]
# [0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 0 0 0]
# [0 0 1 0 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 1 1 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0]
# [0 1 1 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 1 0 0]
# [0 0 0 0 0 0 0 0 1 1 0 1 0 1 0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 1 0]
# [0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 1 1 0 0 0 0 1 1 0 1 0 0 0 1 0 0 0]
# [0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 1 0 0]
# [0 0 0 0 0 0 1 0 0 0 0 0 0 1 1 0 1 0 0 1 0 1 0 1 0 0 0 0 0 0 1 0 0 1 0 0]
# [0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 1 0]
# [0 0 0 1 0 1 0 0 0 1 0 1 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1]
# [0 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0]
# [0 0 0 0 0 0 0 1 0 0 1 0 1 0 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 1 0 1 0 0 0 0]
# [1 0 0 0 0 1 1 1 0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1]
# [0 1 0 0 0 1 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 1 0 0]
# [1 1 1 0 1 1 0 0 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 1 1]
# [1 0 0 0 0 0 0 1 1 0 1 0 0 0 0 0 0 1 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 1 0]
# [0 0 0 0 0 0 0 1 0 0 0 1 1 0 1 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 1 1 0 0 0]
# [0 0 1 0 1 0 0 0 0 0 0 0 1 0 1 0 0 0 1 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 1]
# [1 1 0 0 0 0 0 0 1 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 0 1 0 0]
# [0 0 0 0 0 0 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 0 0]
# [0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 1 0 0 1 1 0 0 0 0 0 1 0 0 1 1 1 0 0 0 0 0]
# [0 0 0 1 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 0 0 0 0 1 1 0 0 1 0 1 0 1 0 0 0 0]
# [1 0 0 0 1 0 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 1 0]
# [0 0 0 0 1 0 1 0 0 1 0 1 0 0 1 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 0 0 0 0 1]
# [0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 1 1 0 0 1 1 0 0 0 0 0 0 0 1 0 0 1 0 0]
# [0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 1 0 1 0 1 0 1 0 0 0 1 0 0 0]
# [0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 1 1 0 0 1 0 0 0 1 0 1 1 0 0 0 0 0 1 0 0]
# [1 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 1]
# [1 0 0 1 0 0 0 0 1 0 0 1 0 1 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1]
# [1 1 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 1 0 0 0]
# [1 1 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 1 1 0 0]
# [0 1 0 0 0 0 0 0 1 1 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 1 0 0 0 1 1 0 0 0]