#Maple program checking the YB relation in type G_2 for the elements h_i^Y(lambda) in Definition 3.1 of our paper "Towards generalized cohomology Schubert calculus via formal root polynomials"

with(combinat);

#the hyperbolic formal group law (a=mu_1, b=mu_2)
f:=proc(x,y) (x+y-a*x*y)/(1+b*x*y) end;

#generates all subsets of {1,...,6} in order to expand the product of h_i(lambda) over the 6 positive roots of G_2
p6:=powerset(6);

#implements Y_i^2=a Y_i
sqr:=proc(s)
local s1, i, coef;
    coef:=1;
    s1:=[s[1]];
    for i from 2 to nops(s) do
        if s[i] <> s[i-1] then s1:=[op(s1), s[i]]
        else coef:=coef*a
        fi
    od;
    [coef, s1]
end;

#calculates a term in the product of h_i(lambda) coresponding to a subset "s" in "p6"
#a product "Y_1.Y_2.Y_1.Y_2" is represented as Y[1,2,1,2], with the appropriate coefficient (product of y_i=y[i], and possibly a, b)
term:=proc(s)
local a, b, i, ys, r;
    if s=[] then 1
    else
        ys:=[];
        for i to nops(s) do
            if s[i] mod 2=0 then ys:=[op(ys), 2]
            else ys:=[op(ys), 1]
            fi
        od;
        r:=sqr(ys);
        (-1)^nops(s)*product('y[s[i]]', 'i'=1 .. nops(s))*r[1]*Y[op(r[2])]
    fi
end;

#calculates a term in the reversed product of h_i(lambda) coresponding to a subset "s" in "p6"
revterm:=proc(s)
local a, b, i, ys, r;
    if s=[] then 1
    else
        ys:=[];
        for i to nops(s) do
            if s[i] mod 2=1 then ys:=[op(ys), 2]
            else ys:=[op(ys), 1]
            fi
        od;
        r:=sqr(ys);
        (-1)^nops(s)*
            product('y[7-s[i]]', 'i'=1 .. nops(s))*r[1]*Y[op(r[2])]
    fi
end;

#expresses the variables y_i in terms of y_1=x and y_6=z corresponding to the simple roots alpha, beta 
#we use the reflection order 
#alpha=(1,-1,0), 
#3 alpha+beta=(2,-1,-1), 
#2 alpha+beta=(1,0,-1), 
#3 alpha+2 beta=(1,1,-2), 
#alpha+beta=(0,1,-1), 
#beta=(-1,2,-1)
sbs:=proc(e)
    simplify(subs({y[1]=x, y[6]=z,
    y[2]=simplify(f(x, f(x, f(x,z)))), 
    y[3]=simplify(f(x, f(x, z))), 
    y[4]=simplify(f(x, f(x, f(x, f(z,z))))),
    y[5]=simplify(f(x,z))}, e))
end;

#calculates the left-hand side in the Yang-Baxter relation
ls:=sum('term([op(p6[i])])','i'=1..nops(p6));

#calculates the right-hand side in the Yang-Baxter relation
rs:=sum('revterm([op(p6[i])])','i'=1..nops(p6));

#calculates the difference of the two sides, after applying the corresponding twisted braid relation in the Demazure algebra (cf. Leclerc, arXiv:1505.05097, Example 4.12 (ii))
dif:=expand(subs(Y[2,1,2,1,2,1]=Y[1,2,1,2,1,2]+4*b*Y[1,2,1,2]-4*b*Y[2,1,2,1]
             +3*b^2*Y[1,2]-3*b^2*Y[2,1],ls-rs));

   y[6] y[2] y[1] a Y[2, 1] - y[6] y[5] y[4] y[1] Y[2, 1, 2, 1]
         + y[6] y[3] y[1] a Y[2, 1]
         - y[6] y[5] y[2] y[1] Y[2, 1, 2, 1]
         + y[5] y[3] y[2] a Y[1, 2] - y[6] y[4] y[3] y[1] a^2  Y[2, 1]
         + y[6] y[5] y[1] a Y[2, 1]
         - y[5] y[4] y[3] y[2] Y[1, 2, 1, 2]
         - 4 y[1] y[2] y[3] y[4] y[5] y[6] b Y[1, 2, 1, 2]
         + 4 y[1] y[2] y[3] y[4] y[5] y[6] b Y[2, 1, 2, 1]
         - 3 y[1] y[2] y[3] y[4] y[5] y[6] b^2  Y[1, 2]
         + 3 y[1] y[2] y[3] y[4] y[5] y[6] b^2  Y[2, 1]
         - y[6] y[5] y[3] y[1] a^2  Y[2, 1]
         - y[6] y[4] y[2] y[1] a^2  Y[2, 1] + y[6] y[5] y[3] a Y[2, 1]
         + y[3] y[4] y[5] y[6] Y[1, 2, 1, 2] - y[6] y[3] Y[2, 1]
         - y[1] y[3] y[4] y[5] y[6] a Y[1, 2, 1, 2]
         + y[5] y[4] y[2] a Y[1, 2]
         + y[1] y[4] y[5] y[6] Y[1, 2, 1, 2]
         - y[1] y[2] y[4] y[5] y[6] a Y[1, 2, 1, 2]
         + y[5] y[6] Y[1, 2] - y[1] y[5] y[6] a Y[1, 2]
         + y[1] y[2] y[5] y[6] Y[1, 2, 1, 2] - y[5] y[4] Y[1, 2]
         - y[3] y[5] y[6] a Y[1, 2] + y[1] y[3] y[5] y[6] a^2  Y[1, 2]
         - y[1] y[2] y[3] y[5] y[6] a Y[1, 2, 1, 2]
         + y[1] y[6] Y[1, 2] - y[1] y[2] y[6] a Y[1, 2]
         + y[3] y[6] Y[1, 2] - y[1] y[3] y[6] a Y[1, 2]
         + y[1] y[2] y[3] y[6] Y[1, 2, 1, 2]
         - y[1] y[4] y[6] a Y[1, 2] + y[1] y[2] y[4] y[6] a^2  Y[1, 2]
         - y[3] y[4] y[6] a Y[1, 2] + y[1] y[3] y[4] y[6] a^2  Y[1, 2]
         - y[1] y[2] y[3] y[4] y[6] a Y[1, 2, 1, 2]
         - y[6] y[5] Y[2, 1] + y[1] y[2] Y[1, 2] + y[2] y[3] Y[2, 1]
         + y[6] y[5] y[4] y[3] y[1] a Y[2, 1, 2, 1]
         + y[1] y[4] Y[1, 2] - y[1] y[2] y[4] a Y[1, 2]
         + y[3] y[4] Y[1, 2] - y[1] y[3] y[4] a Y[1, 2]
         + y[1] y[2] y[3] y[4] Y[1, 2, 1, 2] + y[2] y[5] Y[2, 1]
         - y[6] y[5] y[4] y[3] Y[2, 1, 2, 1]
         - y[2] y[3] y[5] a Y[2, 1] + y[4] y[5] Y[2, 1]
         - y[2] y[4] y[5] a Y[2, 1] + y[6] y[4] y[3] a Y[2, 1]
         + y[2] y[3] y[4] y[5] Y[2, 1, 2, 1] - y[3] y[4] Y[2, 1]
         - y[2] y[3] Y[1, 2] - y[1] y[4] Y[2, 1]
         + y[1] y[3] y[4] a Y[2, 1] - y[1] y[2] Y[2, 1]
         - y[5] y[2] Y[1, 2] + y[1] y[2] y[4] a Y[2, 1]
         - y[1] y[2] y[3] y[4] Y[2, 1, 2, 1]
         + y[6] y[5] y[4] y[2] y[1] a Y[2, 1, 2, 1]
         + y[6] y[4] y[1] a Y[2, 1] - y[6] y[1] Y[2, 1]
         + y[1] y[2] y[3] y[4] y[6] a Y[2, 1, 2, 1]
         - y[1] y[2] y[3] y[6] Y[2, 1, 2, 1]
         + y[1] y[2] y[3] y[5] y[6] a Y[2, 1, 2, 1]

#calculating the difference of the two sides in the Yang-Baxter relation after expressing all the y_i in terms of y_1 and y_6 corresponding to the simple roots
#as we can see from the above expression of "dif", only the coefficients of "Y_1.Y_2.Y_1.Y_2", "Y_2.Y_1.Y_2.Y_1", "Y_1.Y_2", and "Y_2.Y_1" need to be considered; so we only calculate these, and we get "0" each time, as expected.

sbs(coeff(dif,Y[1,2,1,2]));

                                  0

sbs(coeff(dif,Y[2,1,2,1]));

                                  0
sbs(coeff(dif,Y[1,2]));

                                  0

sbs(coeff(dif,Y[2,1]));

                                  0