##########################
#####
#####   Sage calculations for my paper "Residual automorphic forms and spherical unitary representations"
#####   Stephen Miller
#####   Rutgers University
#####   April 2, 2012
#####
#####
#####   The code has 2 parts, the first of which is a Lie file (commented out here).  I saved it as its own file
#####   and ran it from the command line using this command: lie < (filename) > relW_for_E8_node_1 
#####   This then is used in the sage code below
#####
#####



##    LiE code to compute the "Relevant" part of the Weyl group:
##
##    setdefault E8
##    n=8
##    parabolicvector=[1,0,0,0,0,0,0,0]
##    size=W_orbit_size(parabolicvector)
##    print(size)
##    wo=W_orbit(parabolicvector)
##    union(mat a, b) = unique(a^b)
##    intersection(mat a, b) = support(X unique(a)+X unique(b)-X union(a, b))
##    inverse(vec w) = loc l = size(w); loc wi = w;for i = 1 to l do wi[i ] = w[l + 1 - i ] od; wi
##    flipped(vec w) = intersection(pos_roots,W_rt_action(-pos_roots, w))
##    otherroots= support(X id(n)-X parabolicvector)
##    for wor row  wo do if n_rows(intersection(otherroots,flipped(W_word(wor))))!=0 then print(wor) fi od
##    for wor row wo do print(inverse(W_word(wor))) od






#### The sage part:

# This is an encoding of my mathematica program into sage
# A lot of the coding is clumsy, since this is the first project I've tried with sage.

# Reads the file generated from LiE above

file=open('relW_for_E8_node_1','r')
G = WeylGroup("E8",prefix="s");

# sets the parameters 2s and j from the paper.
twoess=8;jay=1;

# Everything below here should be unchanged during a run; only the name of the file, the group name, twoess and jay should be changed.


# Processes the information from LiE in a list called "listfromlie"
# NOTE: this should give an error because I didn't know how to tell it to stop reading.

Lf=file.readlines()
listfromlie=[]
counter=1
for dummy in range(10000000):
    str=Lf[counter]
    while str[-2] != ']':
        counter=counter+1
        str=str+Lf[counter]
    listfromlie=listfromlie+[eval(str)]
    counter=counter+1

# Sage Weyl group and root system data

L = G.domain();s = G.simple_reflections();G.cardinality()
alpha=L.simple_roots();omega=L.fundamental_weights();rho = sum(omega); rho
alphacheck=L.simple_coroots();posroots=L.positive_roots();poscoroots=[2*vec/(vec.scalar(vec)) for vec in posroots]

# Computes W_rel

relW=[];
for entry1 in listfromlie:
    tmp=G.unit();
    for entry2 in entry1:
        tmp=tmp*s[entry2]
    relW=relW+[tmp]


# initializes lambda1 and lambda2 from paper

lambda1=twoess*omega[jay]-rho
lambda2=omega[jay]


# Compute the translates of W_rel

allrelWtranslates=[elt.action(lambda1) for elt in relW]
relWtranslates=uniq(allrelWtranslates)


# Computes which satisfy Langlands' condition (alreadybounded), and which do not.

alreadybounded=[]
notalreadybounded=[]
for elt in relWtranslates:
    tmp=[elt.scalar(om) for om in omega]
    if max(tmp) < 0:
        alreadybounded=alreadybounded+[elt]
    else:
        notalreadybounded=notalreadybounded+[elt]

# We don't use this in the paper, but nonwall is number on the Wall of the negative dual Weyl chamber.

nonwall=[]
for elt in notalreadybounded:
    tmp=[elt.scalar(om) for om in omega]
    if max(tmp)>0:
        nonwall=nonwall+[elt]

# Produces part of the table

print(len(alreadybounded),len(notalreadybounded),len(notalreadybounded)-len(nonwall))


# This tells us which relW elements map lambda1 to a given translate.  It indexes W(lambda1,mu) in the terminology of the paper.

directory=[ [] for i in range(len(relWtranslates))]
for i in  range(len(allrelWtranslates)):
    elt=allrelWtranslates[i]
    spot=relWtranslates.index(elt)
    print(i,spot)
    directory[spot]=directory[spot]+[i]

# computes the involution set

flippedcoroots=[[] for i in range(len(relW))]
for i in range(len(relW)):
    elt=relW[i]
    tmp=[(elt.action(al)).scalar(rho) for al in posroots]
    for entry in range(len(tmp)):
        if tmp[entry] < 0:
            flippedcoroots[i]=flippedcoroots[i]+[entry]

# The t's are used for the Cartan variable in showing nonvanishing (it describes "g")

t1=var('t1')
t2=var('t2')
t3=var('t3')
t4=var('t4')
t5=var('t5')
t6=var('t6')
t7=var('t7')
t8=var('t8')
ts=[t1,t2,t3,t4,t5,t6,t7,t8]

# The inner products of the simple coroots with lambda1 and lambda2

lambda1ips=[lambda1.scalar(pcr) for pcr in poscoroots]
lambda2ips=[lambda2.scalar(pcr) for pcr in poscoroots]


# This is the main data produces by the root system part of the calculation.
# For each W_rel translate mu it computes, for each w in W(lambda1,mu) a 3-tuple.
# The first entry of this 3-tuple is (w lambda_1) acting on H, which is a sum of the t's.
# The second and third entries, respectively, are the inner products of lambda_1 and lambda_2 with positive coroots flipped by w

i1=var('i1')
tableofinnerproducts=[[] for i in range(len(directory))]
for k in range(len(directory)):
    for i in directory[k]:
        out=[[relW[i].action(lambda2).scalar(al) for al in alpha]]
        out=[sum(ts[i1-1]*relW[i].action(lambda2).scalar(alpha[i1]) for i1 in G.index_set())]
        outsub1=[]
        outsub2=[]
        for j in flippedcoroots[i]:
            outsub1=outsub1+[lambda1ips[j]]
            outsub2=outsub2+[lambda2ips[j]]
        out=out+[outsub1,outsub2]
        tableofinnerproducts[k]=tableofinnerproducts[k]+[out]


# This makes a table of the common power of epsilon outside the terms for a fixed mu (and verifies that it is indeed common).

powersofepsilon=[[] for i in range(len(directory))]
for k in range(len(directory)):
    for entry in tableofinnerproducts[k]:
        count=0
        for j in entry[1]:
            if j==1:
                count=count-1
            if j==-1:
                count=count+1
        powersofepsilon[k]=powersofepsilon[k]+[count]
    if len(uniq(powersofepsilon[k])) > 1:
        print("length exception")
    powersofepsilon[k]=powersofepsilon[k][0]


# This prints out those powers, for alreadybounded and notalreadybounded

[powersofepsilon[relWtranslates.index(entry)] for entry in alreadybounded]
[powersofepsilon[relWtranslates.index(entry)] for entry in notalreadybounded]


# Here is a long list of variables to declare

c0p1=var('c0p1')
c0p3=var('c0p3')
c0p5=var('c0p5')
c0p7=var('c0p7')
c0p9=var('c0p9')
c0p11=var('c0p11')
c0p13=var('c0p13')
c0p15=var('c0p15')
c0p17=var('c0p17')
c0p19=var('c0p19')
c0p21=var('c0p21')
c0p23=var('c0p23')
c0p25=var('c0p25')
c0s=[c0p1,c0p3,c0p5,c0p7,c0p9,c0p11,c0p13,c0p15,c0p17,c0p19,c0p21,c0p23,c0p25]
c1p1=var('c1p1')
c1p2=var('c1p2')
c1p3=var('c1p3')
c1p4=var('c1p4')
c1p5=var('c1p5')
c1p6=var('c1p6')
c1p7=var('c1p7')
c1p8=var('c1p8')
c1p9=var('c1p9')
c1p10=var('c1p10')
c1p11=var('c1p11')
c1p12=var('c1p12')
c1p13=var('c1p13')
c1p14=var('c1p14')
c1p15=var('c1p15')
c1p16=var('c1p16')
c1p17=var('c1p17')
c1p18=var('c1p18')
c1p19=var('c1p19')
c1p20=var('c1p20')
c1p21=var('c1p21')
c1p22=var('c1p22')
c1p23=var('c1p23')
c1p24=var('c1p24')
c1p25=var('c1p25')
c1s=[0,c1p1,c1p2,c1p3,c1p4,c1p5,c1p6,c1p7,c1p8,c1p9,c1p10,c1p11,c1p12,c1p13,c1p14,c1p15,c1p16,c1p17,c1p18,c1p19,c1p20,c1p21,c1p22,c1p23,c1p24,c1p25]
c2p0=var('c2p0')
c2p1=var('c2p1')
c2p2=var('c2p2')
c2p3=var('c2p3')
c2p4=var('c2p4')
c2p5=var('c2p5')
c2p6=var('c2p6')
c2p7=var('c2p7')
c2p8=var('c2p8')
c2p9=var('c2p9')
c2p10=var('c2p10')
c2p11=var('c2p11')
c2p12=var('c2p12')
c2p13=var('c2p13')
c2p14=var('c2p14')
c2p15=var('c2p15')
c2p16=var('c2p16')
c2p17=var('c2p17')
c2p18=var('c2p18')
c2p19=var('c2p19')
c2p20=var('c2p20')
c2p21=var('c2p21')
c2p22=var('c2p22')
c2p23=var('c2p23')
c2p24=var('c2p24')
c2p25=var('c2p25')
c2s=[c2p0,c2p1,c2p2,c2p3,c2p4,c2p5,c2p6,c2p7,c2p8,c2p9,c2p10,c2p11,c2p12,c2p13,c2p14,c2p15,c2p16,c2p17,c2p18,c2p19,c2p20,c2p21,c2p22,c2p23,c2p24,c2p25]
c3p0=var('c3p0')
c3p1=var('c3p1')
c3p2=var('c3p2')
c3p3=var('c3p3')
c3p4=var('c3p4')
c3p5=var('c3p5')
c3p6=var('c3p6')
c3p7=var('c3p7')
c3p8=var('c3p8')
c3p9=var('c3p9')
c3p10=var('c3p10')
c3p11=var('c3p11')
c3p12=var('c3p12')
c3p13=var('c3p13')
c3p14=var('c3p14')
c3p15=var('c3p15')
c3p16=var('c3p16')
c3p17=var('c3p17')
c3p18=var('c3p18')
c3p19=var('c3p19')
c3p20=var('c3p20')
c3p21=var('c3p21')
c3p22=var('c3p22')
c3p23=var('c3p23')
c3p24=var('c3p24')
c3p25=var('c3p25')
c3s=[c3p0,c3p1,c3p2,c3p3,c3p4,c3p5,c3p6,c3p7,c3p8,c3p9,c3p10,c3p11,c3p12,c3p13,c3p14,c3p15,c3p16,c3p17,c3p18,c3p19,c3p20,c3p21,c3p22,c3p23,c3p24,c3p25]
c4p0=var('c4p0')
c4p1=var('c4p1')
c4p2=var('c4p2')
c4p3=var('c4p3')
c4p4=var('c4p4')
c4p5=var('c4p5')
c4p6=var('c4p6')
c4p7=var('c4p7')
c4p8=var('c4p8')
c4p9=var('c4p9')
c4p10=var('c4p10')
c4p11=var('c4p11')
c4p12=var('c4p12')
c4p13=var('c4p13')
c4p14=var('c4p14')
c4p15=var('c4p15')
c4p16=var('c4p16')
c4p17=var('c4p17')
c4p18=var('c4p18')
c4p19=var('c4p19')
c4p20=var('c4p20')
c4p21=var('c4p21')
c4p22=var('c4p22')
c4p23=var('c4p23')
c4p24=var('c4p24')
c4p25=var('c4p25')
c4s=[c4p0,c4p1,c4p2,c4p3,c4p4,c4p5,c4p6,c4p7,c4p8,c4p9,c4p10,c4p11,c4p12,c4p13,c4p14,c4p15,c4p16,c4p17,c4p18,c4p19,c4p20,c4p21,c4p22,c4p23,c4p24,c4p25]
c5p0=var('c5p0')
c5p1=var('c5p1')
c5p2=var('c5p2')
c5p3=var('c5p3')
c5p4=var('c5p4')
c5p5=var('c5p5')
c5p6=var('c5p6')
c5p7=var('c5p7')
c5p8=var('c5p8')
c5p9=var('c5p9')
c5p10=var('c5p10')
c5p11=var('c5p11')
c5p12=var('c5p12')
c5p13=var('c5p13')
c5p14=var('c5p14')
c5p15=var('c5p15')
c5p16=var('c5p16')
c5p17=var('c5p17')
c5p18=var('c5p18')
c5p19=var('c5p19')
c5p20=var('c5p20')
c5p21=var('c5p21')
c5p22=var('c5p22')
c5p23=var('c5p23')
c5p24=var('c5p24')
c5p25=var('c5p25')
c5s=[c5p0,c5p1,c5p2,c5p3,c5p4,c5p5,c5p6,c5p7,c5p8,c5p9,c5p10,c5p11,c5p12,c5p13,c5p14,c5p15,c5p16,c5p17,c5p18,c5p19,c5p20,c5p21,c5p22,c5p23,c5p24,c5p25]
c6p0=var('c6p0')
c6p1=var('c6p1')
c6p2=var('c6p2')
c6p3=var('c6p3')
c6p4=var('c6p4')
c6p5=var('c6p5')
c6p6=var('c6p6')
c6p7=var('c6p7')
c6p8=var('c6p8')
c6p9=var('c6p9')
c6p10=var('c6p10')
c6p11=var('c6p11')
c6p12=var('c6p12')
c6p13=var('c6p13')
c6p14=var('c6p14')
c6p15=var('c6p15')
c6p16=var('c6p16')
c6p17=var('c6p17')
c6p18=var('c6p18')
c6p19=var('c6p19')
c6p20=var('c6p20')
c6p21=var('c6p21')
c6p22=var('c6p22')
c6p23=var('c6p23')
c6p24=var('c6p24')
c6p25=var('c6p25')
c6s=[c6p0,c6p1,c6p2,c6p3,c6p4,c6p5,c6p6,c6p7,c6p8,c6p9,c6p10,c6p11,c6p12,c6p13,c6p14,c6p15,c6p16,c6p17,c6p18,c6p19,c6p20,c6p21,c6p22,c6p23,c6p24,c6p25]
c7p0=var('c7p0')
c7p1=var('c7p1')
c7p2=var('c7p2')
c7p3=var('c7p3')
c7p4=var('c7p4')
c7p5=var('c7p5')
c7p6=var('c7p6')
c7p7=var('c7p7')
c7p8=var('c7p8')
c7p9=var('c7p9')
c7p10=var('c7p10')
c7p11=var('c7p11')
c7p12=var('c7p12')
c7p13=var('c7p13')
c7p14=var('c7p14')
c7p15=var('c7p15')
c7p16=var('c7p16')
c7p17=var('c7p17')
c7p18=var('c7p18')
c7p19=var('c7p19')
c7p20=var('c7p20')
c7p21=var('c7p21')
c7p22=var('c7p22')
c7p23=var('c7p23')
c7p24=var('c7p24')
c7p25=var('c7p25')
c7s=[c7p0,c7p1,c7p2,c7p3,c7p4,c7p5,c7p6,c7p7,c7p8,c7p9,c7p10,c7p11,c7p12,c7p13,c7p14,c7p15,c7p16,c7p17,c7p18,c7p19,c7p20,c7p21,c7p22,c7p23,c7p24,c7p25]
c8p0=var('c8p0')
c8p1=var('c8p1')
c8p2=var('c8p2')
c8p3=var('c8p3')
c8p4=var('c8p4')
c8p5=var('c8p5')
c8p6=var('c8p6')
c8p7=var('c8p7')
c8p8=var('c8p8')
c8p9=var('c8p9')
c8p10=var('c8p10')
c8p11=var('c8p11')
c8p12=var('c8p12')
c8p13=var('c8p13')
c8p14=var('c8p14')
c8p15=var('c8p15')
c8p16=var('c8p16')
c8p17=var('c8p17')
c8p18=var('c8p18')
c8p19=var('c8p19')
c8p20=var('c8p20')
c8p21=var('c8p21')
c8p22=var('c8p22')
c8p23=var('c8p23')
c8p24=var('c8p24')
c8p25=var('c8p25')
c8s=[c8p0,c8p1,c8p2,c8p3,c8p4,c8p5,c8p6,c8p7,c8p8,c8p9,c8p10,c8p11,c8p12,c8p13,c8p14,c8p15,c8p16,c8p17,c8p18,c8p19,c8p20,c8p21,c8p22,c8p23,c8p24,c8p25]
c9p0=var('c9p0')
c9p1=var('c9p1')
c9p2=var('c9p2')
c9p3=var('c9p3')
c9p4=var('c9p4')
c9p5=var('c9p5')
c9p6=var('c9p6')
c9p7=var('c9p7')
c9p8=var('c9p8')
c9p9=var('c9p9')
c9p10=var('c9p10')
c9p11=var('c9p11')
c9p12=var('c9p12')
c9p13=var('c9p13')
c9p14=var('c9p14')
c9p15=var('c9p15')
c9p16=var('c9p16')
c9p17=var('c9p17')
c9p18=var('c9p18')
c9p19=var('c9p19')
c9p20=var('c9p20')
c9p21=var('c9p21')
c9p22=var('c9p22')
c9p23=var('c9p23')
c9p24=var('c9p24')
c9p25=var('c9p25')
c9s=[c9p0,c9p1,c9p2,c9p3,c9p4,c9p5,c9p6,c9p7,c9p8,c9p9,c9p10,c9p11,c9p12,c9p13,c9p14,c9p15,c9p16,c9p17,c9p18,c9p19,c9p20,c9p21,c9p22,c9p23,c9p24,c9p25]
c10p0=var('c10p0')
c10p1=var('c10p1')
c10p2=var('c10p2')
c10p3=var('c10p3')
c10p4=var('c10p4')
c10p5=var('c10p5')
c10p6=var('c10p6')
c10p7=var('c10p7')
c10p8=var('c10p8')
c10p9=var('c10p9')
c10p10=var('c10p10')
c10p11=var('c10p11')
c10p12=var('c10p12')
c10p13=var('c10p13')
c10p14=var('c10p14')
c10p15=var('c10p15')
c10p16=var('c10p16')
c10p17=var('c10p17')
c10p18=var('c10p18')
c10p19=var('c10p19')
c10p20=var('c10p20')
c10p21=var('c10p21')
c10p22=var('c10p22')
c10p23=var('c10p23')
c10p24=var('c10p24')
c10p25=var('c10p25')
c10s=[c10p0,c10p1,c10p2,c10p3,c10p4,c10p5,c10p6,c10p7,c10p8,c10p9,c10p10,c10p11,c10p12,c10p13,c10p14,c10p15,c10p16,c10p17,c10p18,c10p19,c10p20,c10p21,c10p22,c10p23,c10p24,c10p25]
c11p0=var('c11p0')
c11p1=var('c11p1')
c11p2=var('c11p2')
c11p3=var('c11p3')
c11p4=var('c11p4')
c11p5=var('c11p5')
c11p6=var('c11p6')
c11p7=var('c11p7')
c11p8=var('c11p8')
c11p9=var('c11p9')
c11p10=var('c11p10')
c11p11=var('c11p11')
c11p12=var('c11p12')
c11p13=var('c11p13')
c11p14=var('c11p14')
c11p15=var('c11p15')
c11p16=var('c11p16')
c11p17=var('c11p17')
c11p18=var('c11p18')
c11p19=var('c11p19')
c11p20=var('c11p20')
c11p21=var('c11p21')
c11p22=var('c11p22')
c11p23=var('c11p23')
c11p24=var('c11p24')
c11p25=var('c11p25')
c11s=[c11p0,c11p1,c11p2,c11p3,c11p4,c11p5,c11p6,c11p7,c11p8,c11p9,c11p10,c11p11,c11p12,c11p13,c11p14,c11p15,c11p16,c11p17,c11p18,c11p19,c11p20,c11p21,c11p22,c11p23,c11p24,c11p25]
c12p0=var('c12p0')
c12p1=var('c12p1')
c12p2=var('c12p2')
c12p3=var('c12p3')
c12p4=var('c12p4')
c12p5=var('c12p5')
c12p6=var('c12p6')
c12p7=var('c12p7')
c12p8=var('c12p8')
c12p9=var('c12p9')
c12p10=var('c12p10')
c12p11=var('c12p11')
c12p12=var('c12p12')
c12p13=var('c12p13')
c12p14=var('c12p14')
c12p15=var('c12p15')
c12p16=var('c12p16')
c12p17=var('c12p17')
c12p18=var('c12p18')
c12p19=var('c12p19')
c12p20=var('c12p20')
c12p21=var('c12p21')
c12p22=var('c12p22')
c12p23=var('c12p23')
c12p24=var('c12p24')
c12p25=var('c12p25')
c12s=[c12p0,c12p1,c12p2,c12p3,c12p4,c12p5,c12p6,c12p7,c12p8,c12p9,c12p10,c12p11,c12p12,c12p13,c12p14,c12p15,c12p16,c12p17,c12p18,c12p19,c12p20,c12p21,c12p22,c12p23,c12p24,c12p25]
c13p0=var('c13p0')
c13p1=var('c13p1')
c13p2=var('c13p2')
c13p3=var('c13p3')
c13p4=var('c13p4')
c13p5=var('c13p5')
c13p6=var('c13p6')
c13p7=var('c13p7')
c13p8=var('c13p8')
c13p9=var('c13p9')
c13p10=var('c13p10')
c13p11=var('c13p11')
c13p12=var('c13p12')
c13p13=var('c13p13')
c13p14=var('c13p14')
c13p15=var('c13p15')
c13p16=var('c13p16')
c13p17=var('c13p17')
c13p18=var('c13p18')
c13p19=var('c13p19')
c13p20=var('c13p20')
c13p21=var('c13p21')
c13p22=var('c13p22')
c13p23=var('c13p23')
c13p24=var('c13p24')
c13p25=var('c13p25')
c13s=[c13p0,c13p1,c13p2,c13p3,c13p4,c13p5,c13p6,c13p7,c13p8,c13p9,c13p10,c13p11,c13p12,c13p13,c13p14,c13p15,c13p16,c13p17,c13p18,c13p19,c13p20,c13p21,c13p22,c13p23,c13p24,c13p25]
c14p0=var('c14p0')
c14p1=var('c14p1')
c14p2=var('c14p2')
c14p3=var('c14p3')
c14p4=var('c14p4')
c14p5=var('c14p5')
c14p6=var('c14p6')
c14p7=var('c14p7')
c14p8=var('c14p8')
c14p9=var('c14p9')
c14p10=var('c14p10')
c14p11=var('c14p11')
c14p12=var('c14p12')
c14p13=var('c14p13')
c14p14=var('c14p14')
c14p15=var('c14p15')
c14p16=var('c14p16')
c14p17=var('c14p17')
c14p18=var('c14p18')
c14p19=var('c14p19')
c14p20=var('c14p20')
c14p21=var('c14p21')
c14p22=var('c14p22')
c14p23=var('c14p23')
c14p24=var('c14p24')
c14p25=var('c14p25')
c14s=[c14p0,c14p1,c14p2,c14p3,c14p4,c14p5,c14p6,c14p7,c14p8,c14p9,c14p10,c14p11,c14p12,c14p13,c14p14,c14p15,c14p16,c14p17,c14p18,c14p19,c14p20,c14p21,c14p22,c14p23,c14p24,c14p25]
c15p0=var('c15p0')
c15p1=var('c15p1')
c15p2=var('c15p2')
c15p3=var('c15p3')
c15p4=var('c15p4')
c15p5=var('c15p5')
c15p6=var('c15p6')
c15p7=var('c15p7')
c15p8=var('c15p8')
c15p9=var('c15p9')
c15p10=var('c15p10')
c15p11=var('c15p11')
c15p12=var('c15p12')
c15p13=var('c15p13')
c15p14=var('c15p14')
c15p15=var('c15p15')
c15p16=var('c15p16')
c15p17=var('c15p17')
c15p18=var('c15p18')
c15p19=var('c15p19')
c15p20=var('c15p20')
c15p21=var('c15p21')
c15p22=var('c15p22')
c15p23=var('c15p23')
c15p24=var('c15p24')
c15p25=var('c15p25')
c15s=[c15p0,c15p1,c15p2,c15p3,c15p4,c15p5,c15p6,c15p7,c15p8,c15p9,c15p10,c15p11,c15p12,c15p13,c15p14,c15p15,c15p16,c15p17,c15p18,c15p19,c15p20,c15p21,c15p22,c15p23,c15p24,c15p25]
c16p0=var('c16p0')
c16p1=var('c16p1')
c16p2=var('c16p2')
c16p3=var('c16p3')
c16p4=var('c16p4')
c16p5=var('c16p5')
c16p6=var('c16p6')
c16p7=var('c16p7')
c16p8=var('c16p8')
c16p9=var('c16p9')
c16p10=var('c16p10')
c16p11=var('c16p11')
c16p12=var('c16p12')
c16p13=var('c16p13')
c16p14=var('c16p14')
c16p15=var('c16p15')
c16p16=var('c16p16')
c16p17=var('c16p17')
c16p18=var('c16p18')
c16p19=var('c16p19')
c16p20=var('c16p20')
c16p21=var('c16p21')
c16p22=var('c16p22')
c16p23=var('c16p23')
c16p24=var('c16p24')
c16p25=var('c16p25')
c16s=[c16p0,c16p1,c16p2,c16p3,c16p4,c16p5,c16p6,c16p7,c16p8,c16p9,c16p10,c16p11,c16p12,c16p13,c16p14,c16p15,c16p16,c16p17,c16p18,c16p19,c16p20,c16p21,c16p22,c16p23,c16p24,c16p25]
c17p0=var('c17p0')
c17p1=var('c17p1')
c17p2=var('c17p2')
c17p3=var('c17p3')
c17p4=var('c17p4')
c17p5=var('c17p5')
c17p6=var('c17p6')
c17p7=var('c17p7')
c17p8=var('c17p8')
c17p9=var('c17p9')
c17p10=var('c17p10')
c17p11=var('c17p11')
c17p12=var('c17p12')
c17p13=var('c17p13')
c17p14=var('c17p14')
c17p15=var('c17p15')
c17p16=var('c17p16')
c17p17=var('c17p17')
c17p18=var('c17p18')
c17p19=var('c17p19')
c17p20=var('c17p20')
c17p21=var('c17p21')
c17p22=var('c17p22')
c17p23=var('c17p23')
c17p24=var('c17p24')
c17p25=var('c17p25')
c17s=[c17p0,c17p1,c17p2,c17p3,c17p4,c17p5,c17p6,c17p7,c17p8,c17p9,c17p10,c17p11,c17p12,c17p13,c17p14,c17p15,c17p16,c17p17,c17p18,c17p19,c17p20,c17p21,c17p22,c17p23,c17p24,c17p25]
c18p0=var('c18p0')
c18p1=var('c18p1')
c18p2=var('c18p2')
c18p3=var('c18p3')
c18p4=var('c18p4')
c18p5=var('c18p5')
c18p6=var('c18p6')
c18p7=var('c18p7')
c18p8=var('c18p8')
c18p9=var('c18p9')
c18p10=var('c18p10')
c18p11=var('c18p11')
c18p12=var('c18p12')
c18p13=var('c18p13')
c18p14=var('c18p14')
c18p15=var('c18p15')
c18p16=var('c18p16')
c18p17=var('c18p17')
c18p18=var('c18p18')
c18p19=var('c18p19')
c18p20=var('c18p20')
c18p21=var('c18p21')
c18p22=var('c18p22')
c18p23=var('c18p23')
c18p24=var('c18p24')
c18p25=var('c18p25')
c18s=[c18p0,c18p1,c18p2,c18p3,c18p4,c18p5,c18p6,c18p7,c18p8,c18p9,c18p10,c18p11,c18p12,c18p13,c18p14,c18p15,c18p16,c18p17,c18p18,c18p19,c18p20,c18p21,c18p22,c18p23,c18p24,c18p25]
c19p0=var('c19p0')
c19p1=var('c19p1')
c19p2=var('c19p2')
c19p3=var('c19p3')
c19p4=var('c19p4')
c19p5=var('c19p5')
c19p6=var('c19p6')
c19p7=var('c19p7')
c19p8=var('c19p8')
c19p9=var('c19p9')
c19p10=var('c19p10')
c19p11=var('c19p11')
c19p12=var('c19p12')
c19p13=var('c19p13')
c19p14=var('c19p14')
c19p15=var('c19p15')
c19p16=var('c19p16')
c19p17=var('c19p17')
c19p18=var('c19p18')
c19p19=var('c19p19')
c19p20=var('c19p20')
c19p21=var('c19p21')
c19p22=var('c19p22')
c19p23=var('c19p23')
c19p24=var('c19p24')
c19p25=var('c19p25')
c19s=[c19p0,c19p1,c19p2,c19p3,c19p4,c19p5,c19p6,c19p7,c19p8,c19p9,c19p10,c19p11,c19p12,c19p13,c19p14,c19p15,c19p16,c19p17,c19p18,c19p19,c19p20,c19p21,c19p22,c19p23,c19p24,c19p25]
c20p0=var('c20p0')
c20p1=var('c20p1')
c20p2=var('c20p2')
c20p3=var('c20p3')
c20p4=var('c20p4')
c20p5=var('c20p5')
c20p6=var('c20p6')
c20p7=var('c20p7')
c20p8=var('c20p8')
c20p9=var('c20p9')
c20p10=var('c20p10')
c20p11=var('c20p11')
c20p12=var('c20p12')
c20p13=var('c20p13')
c20p14=var('c20p14')
c20p15=var('c20p15')
c20p16=var('c20p16')
c20p17=var('c20p17')
c20p18=var('c20p18')
c20p19=var('c20p19')
c20p20=var('c20p20')
c20p21=var('c20p21')
c20p22=var('c20p22')
c20p23=var('c20p23')
c20p24=var('c20p24')
c20p25=var('c20p25')
c20s=[c20p0,c20p1,c20p2,c20p3,c20p4,c20p5,c20p6,c20p7,c20p8,c20p9,c20p10,c20p11,c20p12,c20p13,c20p14,c20p15,c20p16,c20p17,c20p18,c20p19,c20p20,c20p21,c20p22,c20p23,c20p24,c20p25]
c21p0=var('c21p0')
c21p1=var('c21p1')
c21p2=var('c21p2')
c21p3=var('c21p3')
c21p4=var('c21p4')
c21p5=var('c21p5')
c21p6=var('c21p6')
c21p7=var('c21p7')
c21p8=var('c21p8')
c21p9=var('c21p9')
c21p10=var('c21p10')
c21p11=var('c21p11')
c21p12=var('c21p12')
c21p13=var('c21p13')
c21p14=var('c21p14')
c21p15=var('c21p15')
c21p16=var('c21p16')
c21p17=var('c21p17')
c21p18=var('c21p18')
c21p19=var('c21p19')
c21p20=var('c21p20')
c21p21=var('c21p21')
c21p22=var('c21p22')
c21p23=var('c21p23')
c21p24=var('c21p24')
c21p25=var('c21p25')
c21s=[c21p0,c21p1,c21p2,c21p3,c21p4,c21p5,c21p6,c21p7,c21p8,c21p9,c21p10,c21p11,c21p12,c21p13,c21p14,c21p15,c21p16,c21p17,c21p18,c21p19,c21p20,c21p21,c21p22,c21p23,c21p24,c21p25]
c22p0=var('c22p0')
c22p1=var('c22p1')
c22p2=var('c22p2')
c22p3=var('c22p3')
c22p4=var('c22p4')
c22p5=var('c22p5')
c22p6=var('c22p6')
c22p7=var('c22p7')
c22p8=var('c22p8')
c22p9=var('c22p9')
c22p10=var('c22p10')
c22p11=var('c22p11')
c22p12=var('c22p12')
c22p13=var('c22p13')
c22p14=var('c22p14')
c22p15=var('c22p15')
c22p16=var('c22p16')
c22p17=var('c22p17')
c22p18=var('c22p18')
c22p19=var('c22p19')
c22p20=var('c22p20')
c22p21=var('c22p21')
c22p22=var('c22p22')
c22p23=var('c22p23')
c22p24=var('c22p24')
c22p25=var('c22p25')
c22s=[c22p0,c22p1,c22p2,c22p3,c22p4,c22p5,c22p6,c22p7,c22p8,c22p9,c22p10,c22p11,c22p12,c22p13,c22p14,c22p15,c22p16,c22p17,c22p18,c22p19,c22p20,c22p21,c22p22,c22p23,c22p24,c22p25]
c23p0=var('c23p0')
c23p1=var('c23p1')
c23p2=var('c23p2')
c23p3=var('c23p3')
c23p4=var('c23p4')
c23p5=var('c23p5')
c23p6=var('c23p6')
c23p7=var('c23p7')
c23p8=var('c23p8')
c23p9=var('c23p9')
c23p10=var('c23p10')
c23p11=var('c23p11')
c23p12=var('c23p12')
c23p13=var('c23p13')
c23p14=var('c23p14')
c23p15=var('c23p15')
c23p16=var('c23p16')
c23p17=var('c23p17')
c23p18=var('c23p18')
c23p19=var('c23p19')
c23p20=var('c23p20')
c23p21=var('c23p21')
c23p22=var('c23p22')
c23p23=var('c23p23')
c23p24=var('c23p24')
c23p25=var('c23p25')
c23s=[c23p0,c23p1,c23p2,c23p3,c23p4,c23p5,c23p6,c23p7,c23p8,c23p9,c23p10,c23p11,c23p12,c23p13,c23p14,c23p15,c23p16,c23p17,c23p18,c23p19,c23p20,c23p21,c23p22,c23p23,c23p24,c23p25]
c24p0=var('c24p0')
c24p1=var('c24p1')
c24p2=var('c24p2')
c24p3=var('c24p3')
c24p4=var('c24p4')
c24p5=var('c24p5')
c24p6=var('c24p6')
c24p7=var('c24p7')
c24p8=var('c24p8')
c24p9=var('c24p9')
c24p10=var('c24p10')
c24p11=var('c24p11')
c24p12=var('c24p12')
c24p13=var('c24p13')
c24p14=var('c24p14')
c24p15=var('c24p15')
c24p16=var('c24p16')
c24p17=var('c24p17')
c24p18=var('c24p18')
c24p19=var('c24p19')
c24p20=var('c24p20')
c24p21=var('c24p21')
c24p22=var('c24p22')
c24p23=var('c24p23')
c24p24=var('c24p24')
c24p25=var('c24p25')
c24s=[c24p0,c24p1,c24p2,c24p3,c24p4,c24p5,c24p6,c24p7,c24p8,c24p9,c24p10,c24p11,c24p12,c24p13,c24p14,c24p15,c24p16,c24p17,c24p18,c24p19,c24p20,c24p21,c24p22,c24p23,c24p24,c24p25]
c25p0=var('c25p0')
c25p1=var('c25p1')
c25p2=var('c25p2')
c25p3=var('c25p3')
c25p4=var('c25p4')
c25p5=var('c25p5')
c25p6=var('c25p6')
c25p7=var('c25p7')
c25p8=var('c25p8')
c25p9=var('c25p9')
c25p10=var('c25p10')
c25p11=var('c25p11')
c25p12=var('c25p12')
c25p13=var('c25p13')
c25p14=var('c25p14')
c25p15=var('c25p15')
c25p16=var('c25p16')
c25p17=var('c25p17')
c25p18=var('c25p18')
c25p19=var('c25p19')
c25p20=var('c25p20')
c25p21=var('c25p21')
c25p22=var('c25p22')
c25p23=var('c25p23')
c25p24=var('c25p24')
c25p25=var('c25p25')
c25s=[c25p0,c25p1,c25p2,c25p3,c25p4,c25p5,c25p6,c25p7,c25p8,c25p9,c25p10,c25p11,c25p12,c25p13,c25p14,c25p15,c25p16,c25p17,c25p18,c25p19,c25p20,c25p21,c25p22,c25p23,c25p24,c25p25]
c26p0=var('c26p0')
c26p1=var('c26p1')
c26p2=var('c26p2')
c26p3=var('c26p3')
c26p4=var('c26p4')
c26p5=var('c26p5')
c26p6=var('c26p6')
c26p7=var('c26p7')
c26p8=var('c26p8')
c26p9=var('c26p9')
c26p10=var('c26p10')
c26p11=var('c26p11')
c26p12=var('c26p12')
c26p13=var('c26p13')
c26p14=var('c26p14')
c26p15=var('c26p15')
c26p16=var('c26p16')
c26p17=var('c26p17')
c26p18=var('c26p18')
c26p19=var('c26p19')
c26p20=var('c26p20')
c26p21=var('c26p21')
c26p22=var('c26p22')
c26p23=var('c26p23')
c26p24=var('c26p24')
c26p25=var('c26p25')
c26s=[c26p0,c26p1,c26p2,c26p3,c26p4,c26p5,c26p6,c26p7,c26p8,c26p9,c26p10,c26p11,c26p12,c26p13,c26p14,c26p15,c26p16,c26p17,c26p18,c26p19,c26p20,c26p21,c26p22,c26p23,c26p24,c26p25]
c27p0=var('c27p0')
c27p1=var('c27p1')
c27p2=var('c27p2')
c27p3=var('c27p3')
c27p4=var('c27p4')
c27p5=var('c27p5')
c27p6=var('c27p6')
c27p7=var('c27p7')
c27p8=var('c27p8')
c27p9=var('c27p9')
c27p10=var('c27p10')
c27p11=var('c27p11')
c27p12=var('c27p12')
c27p13=var('c27p13')
c27p14=var('c27p14')
c27p15=var('c27p15')
c27p16=var('c27p16')
c27p17=var('c27p17')
c27p18=var('c27p18')
c27p19=var('c27p19')
c27p20=var('c27p20')
c27p21=var('c27p21')
c27p22=var('c27p22')
c27p23=var('c27p23')
c27p24=var('c27p24')
c27p25=var('c27p25')
c27s=[c27p0,c27p1,c27p2,c27p3,c27p4,c27p5,c27p6,c27p7,c27p8,c27p9,c27p10,c27p11,c27p12,c27p13,c27p14,c27p15,c27p16,c27p17,c27p18,c27p19,c27p20,c27p21,c27p22,c27p23,c27p24,c27p25]
c28p0=var('c28p0')
c28p1=var('c28p1')
c28p2=var('c28p2')
c28p3=var('c28p3')
c28p4=var('c28p4')
c28p5=var('c28p5')
c28p6=var('c28p6')
c28p7=var('c28p7')
c28p8=var('c28p8')
c28p9=var('c28p9')
c28p10=var('c28p10')
c28p11=var('c28p11')
c28p12=var('c28p12')
c28p13=var('c28p13')
c28p14=var('c28p14')
c28p15=var('c28p15')
c28p16=var('c28p16')
c28p17=var('c28p17')
c28p18=var('c28p18')
c28p19=var('c28p19')
c28p20=var('c28p20')
c28p21=var('c28p21')
c28p22=var('c28p22')
c28p23=var('c28p23')
c28p24=var('c28p24')
c28p25=var('c28p25')
c28s=[c28p0,c28p1,c28p2,c28p3,c28p4,c28p5,c28p6,c28p7,c28p8,c28p9,c28p10,c28p11,c28p12,c28p13,c28p14,c28p15,c28p16,c28p17,c28p18,c28p19,c28p20,c28p21,c28p22,c28p23,c28p24,c28p25]
c29p0=var('c29p0')
c29p1=var('c29p1')
c29p2=var('c29p2')
c29p3=var('c29p3')
c29p4=var('c29p4')
c29p5=var('c29p5')
c29p6=var('c29p6')
c29p7=var('c29p7')
c29p8=var('c29p8')
c29p9=var('c29p9')
c29p10=var('c29p10')
c29p11=var('c29p11')
c29p12=var('c29p12')
c29p13=var('c29p13')
c29p14=var('c29p14')
c29p15=var('c29p15')
c29p16=var('c29p16')
c29p17=var('c29p17')
c29p18=var('c29p18')
c29p19=var('c29p19')
c29p20=var('c29p20')
c29p21=var('c29p21')
c29p22=var('c29p22')
c29p23=var('c29p23')
c29p24=var('c29p24')
c29p25=var('c29p25')
c29s=[c29p0,c29p1,c29p2,c29p3,c29p4,c29p5,c29p6,c29p7,c29p8,c29p9,c29p10,c29p11,c29p12,c29p13,c29p14,c29p15,c29p16,c29p17,c29p18,c29p19,c29p20,c29p21,c29p22,c29p23,c29p24,c29p25]
c30p0=var('c30p0')
c30p1=var('c30p1')
c30p2=var('c30p2')
c30p3=var('c30p3')
c30p4=var('c30p4')
c30p5=var('c30p5')
c30p6=var('c30p6')
c30p7=var('c30p7')
c30p8=var('c30p8')
c30p9=var('c30p9')
c30p10=var('c30p10')
c30p11=var('c30p11')
c30p12=var('c30p12')
c30p13=var('c30p13')
c30p14=var('c30p14')
c30p15=var('c30p15')
c30p16=var('c30p16')
c30p17=var('c30p17')
c30p18=var('c30p18')
c30p19=var('c30p19')
c30p20=var('c30p20')
c30p21=var('c30p21')
c30p22=var('c30p22')
c30p23=var('c30p23')
c30p24=var('c30p24')
c30p25=var('c30p25')
c30s=[c30p0,c30p1,c30p2,c30p3,c30p4,c30p5,c30p6,c30p7,c30p8,c30p9,c30p10,c30p11,c30p12,c30p13,c30p14,c30p15,c30p16,c30p17,c30p18,c30p19,c30p20,c30p21,c30p22,c30p23,c30p24,c30p25]
cs=[c0s,c1s,c2s,c3s,c4s,c5s,c6s,c7s,c8s,c9s,c10s,c11s,c12s,c13s,c14s,c15s,c16s,c17s,c18s,c19s,c20s,c21s,c22s,c23s,c24s,c25s,c26s,c27s,c28s,c29s,c30s]





# This function takes one of the 3-tuples in tableofinnerproducts and converts it to a pair:
# an outside multiplicative factor times the exponent of an exponential.

def linearfactorandoutsideconstant(tabofinprodentry):
    outsideconstant=1
    linfactor=[0 for i in range(26)]
    linfactor[1]=tabofinprodentry[0]
    for i in range(len(tabofinprodentry[1])):
        if tabofinprodentry[1][i]==0:
            outsideconstant=outsideconstant*-1
            for j in [1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25]:
                linfactor[j]=linfactor[j]+c0s[(j-1)/2]*(tabofinprodentry[2][i])^j
        else:
            if tabofinprodentry[1][i]==-1:
                outsideconstant=outsideconstant*tabofinprodentry[2][i]
                for j in range(1,26):
                    linfactor[j]=linfactor[j]+c1s[j]*(tabofinprodentry[2][i])^j
            else:
                if tabofinprodentry[1][i]==1:
                    outsideconstant=outsideconstant*(-1)/(tabofinprodentry[2][i])
                    for j in range(1,26):
                       linfactor[j]=linfactor[j]-c1s[j]*(-tabofinprodentry[2][i])^j
                else:
                    if tabofinprodentry[1][i] > 1:
                        outsideconstant=outsideconstant*cs[tabofinprodentry[1][i]][0]
                        for j in range(1,26):
                            linfactor[j]=linfactor[j]+cs[tabofinprodentry[1][i]][j]*(tabofinprodentry[2][i])^j
                    else:
                        if tabofinprodentry[1][i] < -1:
                            outsideconstant=outsideconstant/cs[-tabofinprodentry[1][i]][0]
                            for j in range(1,26):
                                linfactor[j]=linfactor[j]-cs[-tabofinprodentry[1][i]][j]*(-tabofinprodentry[2][i])^j
    return [outsideconstant,linfactor]


# More variables to introduce

R.<epsilon>=PowerSeriesRing(SR)
c1=var('c1')
c2=var('c2')
c3=var('c3')
c4=var('c4')
c5=var('c5')
c6=var('c6')
c7=var('c7')
c8=var('c8')
c9=var('c9')
c10=var('c10')
c11=var('c11')
c12=var('c12')
c13=var('c13')
c14=var('c14')
c15=var('c15')
c16=var('c16')
c17=var('c17')
c18=var('c18')
c19=var('c19')
c20=var('c20')

# This is a very slow computation that we can store and reload.  It takes about 8 minutes to do and is common to each calculation.
# It is used for defining the function expseriesfromexponent below
# It is commented out here because I saved it to a file, and we load from it.
#
#expseriestab=exp(R([0,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18,c19,c20],21)).coefficients()
#for j in range(len(expseriestab)):
#    expseriestab[j]=expand(expseriestab[j]*factorial(j))
#save(expseriestab,'expseriesdata')

expseriestab=load('expseriesdata')


# more variables to declare

x1=var('x1')
x2=var('x2')
x3=var('x3')
x4=var('x4')
x5=var('x5')
x6=var('x6')
x7=var('x7')
x8=var('x8')
x9=var('x9')
x10=var('x10')
x11=var('x11')
x12=var('x12')
x13=var('x13')
x14=var('x14')
x15=var('x15')
x16=var('x16')
x17=var('x17')
x18=var('x18')
x19=var('x19')
x20=var('x20')

# Thus function evaluates derivatives of exp(x1*x+x2*x^2+...) at x=0

def expseriesfromexponent(inputvector,power):
    return expseriestab[power].substitute(c1=inputvector[1]).substitute(c2=inputvector[2]).substitute(c3=inputvector[3]).substitute(c4=inputvector[4]).substitute(c5=inputvector[5]).substitute(c6=inputvector[6]).substitute(c7=inputvector[7]).substitute(c8=inputvector[8]).substitute(c9=inputvector[9]).substitute(c10=inputvector[10]).substitute(c11=inputvector[11]).substitute(c12=inputvector[12]).substitute(c13=inputvector[13]).substitute(c14=inputvector[14]).substitute(c15=inputvector[15]).substitute(c16=inputvector[16]).substitute(c17=inputvector[17]).substitute(c18=inputvector[18]).substitute(c19=inputvector[19]).substitute(c20=inputvector[20])

# This command calculates a series coefficient for a particular mu, but does not take into account the overall power of epsilon.
# Thus power is nonnegative here.  We normalize by setting the first outside constant to 1, to simplify the expression.

def preC(muindex,power):
    rawsource=map(linearfactorandoutsideconstant, tableofinnerproducts[muindex])
    output=0
    for entry in rawsource:
        output=output+entry[0]*expseriesfromexponent(entry[1],power)/rawsource[0][0]
    return simplify(expand(output))

# This is C(mu,k,g) from the paper.  g is indexed by the t's.

def C(muindex,kay):
    if kay-powersofepsilon[muindex] >= 0: return preC(muindex,kay-powersofepsilon[muindex])
    if kay-powersofepsilon[muindex] < 0: return 0

# smallest power of epsilon in calculation

emin=min(powersofepsilon)

# To check alreadybounded terms vanish to order < ord

def alreadyboundedchecker(ord):
    return [[C(relWtranslates.index(entry),k) for k in range(emin,ord)] for entry in alreadybounded]


# To check notalreadybounded terms vanish to order < ord

def notalreadyboundedchecker(ord):
    return [[C(relWtranslates.index(entry),k) for k in range(emin,ord)] for entry in notalreadybounded]


#    ---------

# alreadyboundedchecker with printing counter:

def alreadyboundedcheckerprint(ord,start,end):
    output=[]
    for j in range(start,end):
        if Mod(j,100)==0: print(j)
        output=output+[C(j,k) for k in range(emin,ord)]
    return output

# To hunt for a nonvanishing alreadybounded term:

C(relWtranslates.index(alreadybounded[0]),1)