# Construction of function field of prounipotent mixed Tate Galois group of ZZ[1/6] as polynomial algebra on Lyndon words. We replace the Article's tau, upsilon, sigma by the variables a, b, and c, respectively. 

AZ.<a, b, ab, c, aab, abb, bc, ac, aaab, aabb, abbb>=QQ[]

KZ = Frac(AZ)

#construction of words in KZ

#half-weight 2

aa = (1/2)*a*a
ba = a*b - ab
bb = (1/2)*b*b

#half-weight 3

aaa = (1/6)*a*a*a
bbb = (1/6)*b*b*b

aba = a*ab - 2*aab
baa = b*aa - aba - aab
bab = 2*bb*a - b*ba - 2*abb
bba = b*ba - bb*a + abb

#half-weight 4

ca = a*c - ac
cb = b*c - bc

aaaa = (1/24)*a*a*a*a
bbbb = (1/24)*b*b*b*b

aaba = a*aab - 3*aaab
abaa = ab*aa - 2*aaba - 3*aaab
baaa = b*aaa - abaa - aaba - aaab
abab = (1/2)*ab*ab - 2*aabb
abba = abb*a - 2*aabb - abab
baab = b*aab - abab - 2*aabb
baba = b*aba - 2*abba - abab
bbaa = (1/2)*b*baa - (1/2)*baba - (1/2)*baab
babb = b*abb - 3*abbb
bbab = bb*ab -  2*babb - 3*abbb
bbba = (1/3)*b*bba - (1/3)*bbab

# Construction of coordinate ring of Albanese variety (= polylog quotient of unipotent fundamental group of X). 

AX.<fe1,fe2,fe11,fe22,fe12,fe11e1,fe22e2,fe12e3,fe11e1e1,fe22e2e2,fe12e3e3,fe11e1e1e1,fe22e2e2e2,fe12e3e3e3> = KZ[]
AXweights =\
 [1,1,1,1,1,2,2,2,3,3,3,4,4,4]

fe3 = fe1 + fe2

#Construction of F

AX_Y.<Y> = AX[]

AX_xy.<x,y> = AX[]

A1 = -(ab - (1/2)*a*b) * (bb * fe11)

B1 = (ab - (1/2)*a*b) * ( b *fe11e1 - ba*Y *fe11 ) + (ba - (1/2) *b*a) * (ab* fe11) * Y -  (ab*b - a*bb) * ( fe11e1 - (1/2)* fe11* fe1 )

C1 = (ba - (1/2)*b*a) * ( -a*fe11e1 + aa*Y* fe11 )* Y - (aa*b - a*bb) *Y*(fe11e1 - (1/2)*fe11*fe1 )

Q1 = - ( (ab*b - a*bb)*x + ( aa*b - a*bb )*y )*fe11e1e1e1  + ( b *fe11e1 - (ba*y + bb*x)*fe11 ) *(aaaa*y^3 + abbb*x^3 + (aaab + aaba + abaa) * x*y^2 + (abba + abab + aabb)*y*x^2 ) + (-a* fe11e1 + (ab*x + aa*y)*fe11)*(baaa*y^3 + bbbb*x^3 + (baab + baba + bbaa) * x*y^2  + (bbba + bbab + babb) *y*x^2 ) + (ca* y +  cb* x)*(1/c)* ( fe11e1e1 - (b*fe11e1 - (ba*y + bb*x)*fe11)*(aaa*y^2 + abb*x^2 + aab*x*y + aba * x*y) + (-a*fe11e1 + (ab*x + aa*y)* fe11)*(baa*y^2 + bbb*x^2 + bab*x*y + bba*x*y) )

Q1 = AX_xy(Q1)

def a1(i,j):
	return Q1.monomial_coefficient(x^i*y^j)

def snd11(i,k):
	return binomial(i,2*k)*(B1^2-4*A1*C1)^k * (-B1)^(i-2*k)

def s1(i):
	return sum(snd11(i,k) for k in range(0,floor(i/2)+1))

def s21(i):
	return sum(a1(i,j)*(2*A1)^(4-i)*Y^j*s1(i) for j in range(0,5))

Qa1 = sum(s21(i) for i in range(0,5))

def bsnd11(i,k):
	return binomial(i,2*k+1)*(B1^2-4*A1*C1)^k * (-B1)^(i-2*k-1)

def bs11(i):
	return sum(bsnd11(i,k) for k in range(0,floor((i-1)/2)+1))

def bs21(i):
	return sum(a1(i,j)*(2*A1)^(4-i)*Y^j*bs11(i) for j in range(0,5))

Qb1 = sum(bs21(i) for i in range(0,5))

P1 = Qa1^2 - (B1^2-4*A1*C1)*Qb1^2

A2 = -(ab - (1/2)*a*b) * (bb * fe22)

B2 = (ab - (1/2)*a*b) * ( b *fe22e2 - ba*Y *fe22 ) + (ba - (1/2) *b*a) * (ab* fe22) * Y -  (ab*b - a*bb) * ( fe22e2 - (1/2)* fe22* fe2 )

C2 = (ba - (1/2)*b*a) * ( -a*fe22e2 + aa*Y* fe22 )* Y - (aa*b - a*bb) *Y*(fe22e2 - (1/2)*fe22*fe2 )

Q2 = - ( (ab*b - a*bb)*x + ( aa*b - a*bb )*y )*fe22e2e2e2  + ( b *fe22e2 - (ba*y + bb*x)*fe22 ) *(aaaa*y^3 + abbb*x^3 + (aaab + aaba + abaa) * x*y^2 + (abba + abab + aabb)*y*x^2 ) + (-a* fe22e2 + (ab*x + aa*y)*fe22)*(baaa*y^3 + bbbb*x^3 + (baab + baba + bbaa) * x*y^2  + (bbba + bbab + babb) *y*x^2 ) + (ca* y +  cb* x)*(1/c)* ( fe22e2e2 - (b*fe22e2 - (ba*y + bb*x)*fe22)*(aaa*y^2 + abb*x^2 + aab*x*y + aba * x*y) + (-a*fe22e2 + (ab*x + aa*y)* fe22)*(baa*y^2 + bbb*x^2 + bab*x*y + bba*x*y) )

Q2 = AX_xy(Q2)

def a2(i,j):
	return Q2.monomial_coefficient(x^i*y^j)

def snd12(i,k):
	return binomial(i,2*k)*(B2^2-4*A2*C2)^k * (-B2)^(i-2*k)

def s2(i):
	return sum(snd12(i,k) for k in range(0,floor(i/2)+1))

def s22(i):
	return sum(a2(i,j)*(2*A2)^(4-i)*Y^j*s2(i) for j in range(0,5))

Qa2 = sum(s22(i) for i in range(0,5))

def bsnd12(i,k):
	return binomial(i,2*k+1)*(B2^2-4*A2*C2)^k * (-B2)^(i-2*k-1)

def bs12(i):
	return sum(bsnd12(i,k) for k in range(0,floor((i-1)/2)+1))

def bs22(i):
	return sum(a2(i,j)*(2*A2)^(4-i)*Y^j*bs12(i) for j in range(0,5))

Qb2 = sum(bs22(i) for i in range(0,5))

P2 = Qa2^2 - (B2^2-4*A2*C2)*Qb2^2

A3 = -(ab - (1/2)*a*b) * (bb * fe12)

B3 = (ab - (1/2)*a*b) * ( b *fe12e3 - ba*Y *fe12 ) + (ba - (1/2) *b*a) * (ab* fe12) * Y -  (ab*b - a*bb) * ( fe12e3 - (1/2)* fe12* fe3 )

C3 = (ba - (1/2)*b*a) * ( -a*fe12e3 + aa*Y* fe12 )* Y - (aa*b - a*bb) *Y*(fe12e3 - (1/2)*fe12*fe3 )

Q3 = - ( (ab*b - a*bb)*x + ( aa*b - a*bb )*y )*fe12e3e3e3  + ( b *fe12e3 - (ba*y + bb*x)*fe12 ) *(aaaa*y^3 + abbb*x^3 + (aaab + aaba + abaa) * x*y^2 + (abba + abab + aabb)*y*x^2 ) + (-a* fe12e3 + (ab*x + aa*y)*fe12)*(baaa*y^3 + bbbb*x^3 + (baab + baba + bbaa) * x*y^2  + (bbba + bbab + babb) *y*x^2 ) + (ca* y +  cb* x)*(1/c)* ( fe12e3e3 - (b*fe12e3 - (ba*y + bb*x)*fe12)*(aaa*y^2 + abb*x^2 + aab*x*y + aba * x*y) + (-a*fe12e3 + (ab*x + aa*y)* fe12)*(baa*y^2 + bbb*x^2 + bab*x*y + bba*x*y) )

Q3 = AX_xy(Q3)

def a3(i,j):
	return Q3.monomial_coefficient(x^i*y^j)

def snd13(i,k):
	return binomial(i,2*k)*(B3^2-4*A3*C3)^k * (-B3)^(i-2*k)

def s3(i):
	return sum(snd13(i,k) for k in range(0,floor(i/2)+1))

def s23(i):
	return sum(a3(i,j)*(2*A3)^(4-i)*Y^j*s3(i) for j in range(0,5))

Qa3 = sum(s23(i) for i in range(0,5))

def bsnd13(i,k):
	return binomial(i,2*k+1)*(B3^2-4*A3*C3)^k * (-B3)^(i-2*k-1)

def bs13(i):
	return sum(bsnd13(i,k) for k in range(0,floor((i-1)/2)+1))

def bs23(i):
	return sum(a3(i,j)*(2*A3)^(4-i)*Y^j*bs13(i) for j in range(0,5))

Qb3 = sum(bs23(i) for i in range(0,5))

P3 = Qa3^2 - (B3^2-4*A3*C3)*Qb3^2

p1,trash = P1.quo_rem(Y)

p2,trash = P2.quo_rem(Y)

p3,trash = P3.quo_rem(Y)

AX_X_Y.<X> = AX_Y[]

p1X = AX_X_Y(p1(X))

p2YX = p2(Y-X)

Res12 = p1X.resultant(p2YX)

F = Res12.resultant(p3)