# Expand the log of the KZ Drinfeld associator up to a certain degree
# over the Hall bases
# Needs: mzv-1-12.txt (the French tables)
#        bas_comm.2 -- bas_comm.N
#        assolist.2 -- N
# no subtitution for zeta(2)!
# 
# version F: expanding over the basic commutators, not over
# the basic zeta numbers, and substituting numerical values of zetas

Digits:=100;

interface(screenwidth=350);
interface(echo=0);
interface(prettyprint=0);
interface(quiet=true);
writeto(`asso12f`):

with(linalg):
read(`mzv-1-12w.txt`):
read(`assolist.2-8`):
read(`assolist.09`):
read(`assolist.10`):

# Concatenation multiplication:
cm:= proc(p1,p2)
local e1,e2,t1,t2,v1,v2,n1,n2,i1,i2,c1,c2,v,r;
e1:=expand(p1);
e2:=expand(p2);
r:=0;
for i1 from 1 to nops(e1) do
for i2 from 1 to nops(e2) do
t1:=op(i1,e1);
t2:=op(i2,e2);
n1:=nops(t1);
n2:=nops(t2);
v1:=op(n1,t1);
v2:=op(n2,t2);
c1:=coeff(t1,v1);
c2:=coeff(t2,v2);
r:=r+c1*c2*cat(v1,v2);
od;od;
RETURN(r);
end:
# Convert a linear combination of variables
# AB,BA,AAB,ABA,ABB,BAA,BAB,BBA,AAAB... into a vector
digit01:=proc(a)
if a=65 then 0 else 1 fi;
end:
numvar:= proc(a)
local la,i,r;
la:=convert(a,bytes);
r:=0;
for i from 1 to nops(la) do
r:=r*2+digit01(la[i]); od;
RETURN(r);
end:
tovect:=proc(n,p)
local a,k,i,c,L,v;
L:=2^n-2;
for i from 1 to L do v[i]:=0 od;
for i from 1 to nops(p) do
a:=op(i,p);
if nops(a)=1 then c:=1;k:=numvar(a)
else c:=op(1,a);k:=numvar(op(2,a)) fi;
v[k]:=v[k]+c;
od;
RETURN([seq(v[i],i=1..L)]);
end:

coeffm:= proc(p,v)
local t,v1,v2,n,r;
t:=whattype(v);
if t='symbol' then r:=coeff(p,v)
elif t='function' then r:=coeff(p,v)
elif t='`*`' then 
v1:=op(1,v); v2:=product(op(ii,v),ii=2..nops(v)); r:=coeffm(coeff(p,v1),v2)
elif t='`^`' then 
v1:=op(1,v); n:=op(2,v); r:=coeff(p,v1,n)
fi;
RETURN(r);
end;

read(`assolist.12`):

logphi12:=ass12

-1/2*(cm(ass2,ass10)+cm(ass3,ass9)+cm(ass4,ass8)+cm(ass5,ass7)
+cm(ass6,ass6)+cm(ass7,ass5)+cm(ass8,ass4)+cm(ass9,ass3)+cm(ass10,ass2))

+1/3*(
cm(cm(ass2,ass2),ass8)+cm(ass2,cm(ass8,ass2))+cm(ass8,cm(ass2,ass2))

+cm(cm(ass2,ass3),ass7)+cm(cm(ass2,ass7),ass3)+cm(cm(ass3,ass2),ass7)
+cm(cm(ass3,ass7),ass2)+cm(cm(ass7,ass2),ass3)+cm(cm(ass7,ass3),ass2)

+cm(cm(ass2,ass4),ass6)+cm(cm(ass2,ass6),ass4)+cm(cm(ass4,ass2),ass6)
+cm(cm(ass4,ass6),ass2)+cm(cm(ass6,ass2),ass4)+cm(cm(ass6,ass4),ass2)

+cm(ass2,cm(ass5,ass5))+cm(ass5,cm(ass2,ass5))+cm(ass5,cm(ass5,ass2))

+cm(ass3,cm(ass3,ass6))+cm(ass3,cm(ass6,ass3))+cm(ass6,cm(ass3,ass3))

+cm(ass3,cm(ass4,ass5))+cm(ass3,cm(ass5,ass4))+cm(ass4,cm(ass3,ass5))
+cm(ass4,cm(ass5,ass3))+cm(ass5,cm(ass3,ass4))+cm(ass5,cm(ass4,ass3))

+cm(ass4,cm(ass4,ass4))
)

-1/4*(
cm(ass2,cm(ass2,cm(ass2,ass6)))+cm(ass2,cm(ass2,cm(ass6,ass2)))
+cm(ass2,cm(ass6,cm(ass2,ass2)))+cm(ass6,cm(ass2,cm(ass2,ass2)))

+cm(ass2,cm(ass2,cm(ass3,ass5)))+cm(ass2,cm(ass3,cm(ass2,ass5))) 
+cm(ass2,cm(ass3,cm(ass5,ass2)))+cm(ass3,cm(ass2,cm(ass2,ass5)))
+cm(ass3,cm(ass2,cm(ass5,ass2)))+cm(ass3,cm(ass5,cm(ass2,ass2)))
+cm(ass2,cm(ass2,cm(ass5,ass3)))+cm(ass2,cm(ass5,cm(ass2,ass3))) 
+cm(ass2,cm(ass5,cm(ass3,ass2)))+cm(ass5,cm(ass2,cm(ass2,ass3)))
+cm(ass5,cm(ass2,cm(ass3,ass2)))+cm(ass5,cm(ass3,cm(ass2,ass2)))

+cm(ass2,cm(ass2,cm(ass4,ass4)))+cm(ass2,cm(ass4,cm(ass2,ass4)))
+cm(ass2,cm(ass4,cm(ass4,ass2)))+cm(ass4,cm(ass2,cm(ass2,ass4)))
+cm(ass4,cm(ass2,cm(ass4,ass2)))+cm(ass4,cm(ass4,cm(ass2,ass2)))

+cm(ass2,cm(ass3,cm(ass3,ass4)))+cm(ass2,cm(ass3,cm(ass4,ass3)))
+cm(ass2,cm(ass4,cm(ass3,ass3)))+cm(ass3,cm(ass2,cm(ass3,ass4)))
+cm(ass3,cm(ass2,cm(ass4,ass3)))+cm(ass4,cm(ass2,cm(ass3,ass3)))
+cm(ass3,cm(ass3,cm(ass2,ass4)))+cm(ass3,cm(ass4,cm(ass2,ass3)))
+cm(ass4,cm(ass3,cm(ass2,ass3)))+cm(ass3,cm(ass3,cm(ass4,ass2)))
+cm(ass3,cm(ass4,cm(ass3,ass2)))+cm(ass4,cm(ass3,cm(ass3,ass2)))

+cm(ass3,cm(ass3,cm(ass3,ass3)))
)

+1/5*(
cm(ass2,cm(ass2,cm(ass2,cm(ass2,ass4))))
+cm(ass2,cm(ass2,cm(ass2,cm(ass4,ass2))))
+cm(ass2,cm(ass2,cm(ass4,cm(ass2,ass2))))
+cm(ass2,cm(ass4,cm(ass2,cm(ass2,ass2))))
+cm(ass4,cm(ass2,cm(ass2,cm(ass2,ass2))))

+cm(ass2,cm(ass2,cm(ass2,cm(ass3,ass3))))
+cm(ass2,cm(ass2,cm(ass3,cm(ass2,ass3))))
+cm(ass2,cm(ass3,cm(ass2,cm(ass2,ass3))))
+cm(ass3,cm(ass2,cm(ass2,cm(ass2,ass3))))
+cm(ass2,cm(ass2,cm(ass3,cm(ass3,ass2))))
+cm(ass2,cm(ass3,cm(ass2,cm(ass3,ass2))))
+cm(ass3,cm(ass2,cm(ass2,cm(ass3,ass2))))
+cm(ass2,cm(ass3,cm(ass3,cm(ass2,ass2))))
+cm(ass3,cm(ass2,cm(ass3,cm(ass2,ass2))))
+cm(ass3,cm(ass3,cm(ass2,cm(ass2,ass2))))
)

-1/6*cm(ass2,cm(ass2,cm(ass2,cm(ass2,cm(ass2,ass2)))))
:
listz:=[z2^6, z2^3*z3^2, z3^4, z2^2*z3*z5, z2*z5^2, z2^2*z62, 
      z2*z3*z7, z5*z7, z3*z9, z2*z82, z8211, zA2];
len:=nops(listz);

read(`bas_comm.12`):
matr12t:=transpose(matr12):

i:='i';
for i from 1 to len do
log12z||i:=simplify(coeffm(logphi12,listz[i])):
od:

j:='j':
difference:=simplify(logphi12-sum('log12z||j'*listz[j],j=1..len)):
if difference=0 then print('difference=0') fi;
j:='j':

i:='i':
for i from 1 to len do
n:=numer(log12z||i):
d:=denom(log12z||i):
c[i]:=listz[i]/d:
nv:=tovect(12,n):
v[i]:=linsolve(matr12t,nv):
od:

size:=vectdim(matr12);
read(`z_numer.mpl`);
for j from 1 to size do lprint (sum(c[k]*v[k][j],k=1..len),"Y") od;