> # This program computes the cocycle invariant 
> # for the 3-twist spun trefoil tau^3(T(2,3)),
> # for the quandle S_4, with coefficient Z_2.
> restart;
> with(linalg):
> 
Warning, new definition for norm
Warning, new definition for trace
>  # Initial conditions:
>  p :=2: #mod p LX
>  q :=2: #A=Z_q
>  Rel0:=T^2+T+1: #order ideal
>  Rel1:=T^2=-T-1: # ideal solved for T^2
>  Rel2:=T^(-1)=-T-1: #inverse rel from Rel1
> 
> # Quandle operation:
> g:=proc(a, b)
> algsubs( Rel1,  
> expand( T*a + (1-T)*b) ) mod p :  end: 
> 
> # Conversion between integers and quandle elements : 
> f:=array[0..p^2-1]:
>    for jj from 0 to p^2-1 do
>    jj1:=iquo(jj,p):
>    jj2:=irem(jj,p):
>    f[jj]:=jj2+jj1*T:
>    od:
> 
>  finv:=proc(c)
>    (coeff(c, T, 0) mod p) 
>   +(coeff(c, T, 1) mod p)*p  : end:
> 
> # A choice of cocycles, made after experiments:
> assign(c[45]=1):
> assign(c[46]=0):
> assign(c[47]=0):
> assign(c[53]=0):
> assign(c[55]=0):
> assign(c[56]=0):
> assign(c[57]=0):
> assign(c[58]=0):
> assign(c[60]=1):
> 
> # Defining a matrix of cocycle conditions:
> A:=array(sparse,1..(p^8+2*p^4),1..p^6):
>  
>  ind:=0:
>  
>  for x from 0 to p^2-1 do
>  for y from 0 to p^2-1 do
>  for z from 0 to p^2-1 do
>  for s from 0 to p^2-1 do
>  
>  xstry:=coeff(g(f[x],f[y]),T,0)+coeff(g(f[x],f[y]),T,1)*p:
>  xstrz:=coeff(g(f[x],f[z]),T,0)+coeff(g(f[x],f[z]),T,1)*p:
>  xstrs:=coeff(g(f[x],f[s]),T,0)+coeff(g(f[x],f[s]),T,1)*p:
>  ystrz:=coeff(g(f[y],f[z]),T,0)+coeff(g(f[y],f[z]),T,1)*p:
>  ystrs:=coeff(g(f[y],f[s]),T,0)+coeff(g(f[y],f[s]),T,1)*p:
>  zstrs:=coeff(g(f[z],f[s]),T,0)+coeff(g(f[z],f[s]),T,1)*p:
>  
>  ind:=ind+1:
> 
>  if abs(x-y)>0 then
>  if abs(y-z)>0 then
>   A[ind,x*p^4 + y*p^2 + z +1]:= 
>   (1+A[ind,x*p^4 +  y*p^2 + z +1]) mod q:
>  fi:
>  fi:
>  
>  if abs(x-z) > 0 then
>  if abs(z-s) > 0 then
>   A[ind,x*p^4 + z*p^2 +s+1] := 
> (1 + A[ind,x*p^4 + z*p^2 +s+1]) mod q:
>  fi: 
>  fi:
>  
>  if abs(xstrz-ystrz) >0 then
>  if abs(ystrz-s) >0 then
>   A[ind,xstrz*p^4 + ystrz*p^2 +s+1]:= 
> (1+ A[ind,xstrz*p^4 + ystrz*p^2+s+1]) mod q:
>  fi:
>  fi:
>  
>  if abs(xstry-z) >0 then
>  if abs(z-s) >0 then
>    A[ind,xstry*p^4 + z* p^2 + s+1]:=
>  (-1+A[ind,xstry*p^4 + z* p^2 + s+1]) mod q:
>  fi:
>  fi:
>  
>  if abs(x-y) >0 then
>  if abs(y-s) >0 then
>   A[ind,x*p^4+y*p^2+s+1]:=
>   (-1+A[ind,x*p^4+y*p^2+s+1]) mod q:
>  fi:
>  fi:
>  
>  if abs(xstrs-ystrs) > 0 then
>  if abs(ystrs-zstrs) > 0 then
> 
> A[ind,xstrs*p^4+ystrs*p^2+zstrs+1]:=
> (-1+A[ind,xstrs*p^4+ystrs*p^2+zstrs+1]) mod q:
>  fi:
>  fi:
>  od: od: od: od:
>  
>  for w1 from 0 to p^2-1 do
>  for w2 from 0 to p^2-1 do
>  A[p^6 + w1*p^2 +w2 +1, w1*p^4 + w1*p^2 +w2 +1]:=1:
>  A[p^6 + p^4+w1*p^2+w2+1,w1*p^4+w2*p^2+w2+1]:=1:
>  od:
>  od:
>  
> # Solving the cocycle conditions:
>  Sol:=Linsolve(A,vector(p^8+2*p^4,0),'rank',c) mod q : 
>  # 
> # Writing the state-sum expressions:
> # (This was obtained in the paper 
> # ``Computations of quandle cocycles....'')
> for x from 0 to p^2-1 do
> for y from 0 to p^2-1 do
> 
> Theta03[x,y]:= 
>  Sol[ finv( algsubs(Rel2, 
>       expand( T^(-1)*f[y]+(1-T^(-1))*f[x] ) mod p ))*p^4 
>     + (x)*p^2 
>     + (y) +1 ] +
>  Sol[ finv( algsubs(Rel2, 
>       expand( T^(-1)*f[y]+(1-T^(-1))*f[x] ) mod p  ))*p^4     +
> (y)*p^2
>     + finv( algsubs(Rel1,
>       expand( T*f[x]+(1-T)*f[y] ) mod p )) +1 ] :
> ### the above line should be 
> ### expand( T*f[x]-(1-T)*f[y] ) mod p )) +1 ] :
> ### (- instead of + , see the paper ``Computations...''
> ### and this typo appears in some other programs. 
> ### But forturnately, +=- in case p=2, and all files
> ### with this typo had p=2 ! 
> 
> Theta13[x,y]:= 
>  Sol[ finv( algsubs(Rel2, 
>       expand( T^(-1)*f[y]+(1-T^(-1))*f[x] mod p ) ))*p^4 
>     + (x)*p^2
>     + finv( algsubs(Rel2, 
>       expand( T^(-1)*f[y]+(1-T^(-1))*f[x] ) mod p ) )+1 ] +Sol[
> (x)*p^4
>     + (y)*p^2
>     + finv( algsubs(Rel2, 
>       expand( T^(-1)*f[y]+(1-T^(-1))*f[x] ) 
>       mod p ) ) +1 ] :
> 
> od: od: 
> 
> 
> $ Evaluating the state-sum:
> S:=0:
> 
> for x from 0 to p^2-1 do 
> for y from 0 to p^2-1 do 
> 
> E:= - Theta03[x,y] 
>     - Theta03[x,
>       finv(algsubs(Rel2,
>       expand(T^(-1)*f[y]-(1-T^(-1))*f[x])
>       mod p ) ) ]
>     + Theta13[x,y]
>     + Theta13[x, 
>       finv(algsubs(Rel2,
>       expand(T^(-1)*f[y]-(1-T^(-1))*f[x])
>       mod p ) ) ] mod q: 
> 
> S:= S+ t^E :
> 
> od: od:
> 
> print(S);
> 
> 

                               8 + 8 t

> # Finding the cocycles that was used:
> cocy:=0:
> 
> for sub1 from 1 to p^6 do 
> cocy:=cocy + Sol[sub1]*(chi[f[iquo(sub1-1,p^4)],
>            f[iquo(irem(sub1-1, p^4),p^2)], 
>            f[irem(sub1-1,p^2)] ] ) :
> od:
> 
> print(cocy);

  c[31] chi[1, 1 + T, 0] + chi[1, T, 1 + T] + chi[1, T, 0]

         + chi[1, T, 1] + c[31] chi[1, 0, 1 + T]

         + (1 + c[31]) chi[1, 0, 1] + c[40] chi[0, 1 + T, T]

         + (c[40] + c[31]) chi[0, T, 1 + T] + c[40] chi[0, 1 + T, 1]

         + chi[0, 1 + T, 0] + (1 + c[40]) chi[0, T, 0]

         + (1 + c[31]) chi[0, T, 1] + c[31] chi[0, 1, 0]

         + (1 + c[31]) chi[0, 1, T] + chi[T, 1 + T, 0]

         + c[40] chi[T, 1, 1 + T] + chi[T, 1, T] + c[40] chi[T, 1, 0]

         + (1 + c[40]) chi[T, 0, T] + chi[T, 0, 1 + T]

         + c[40] chi[T, 0, 1] + c[31] chi[1, 1 + T, T]

         + chi[1 + T, T, 1 + T] + chi[1 + T, 0, T]

         + chi[1 + T, 0, 1 + T] + chi[1 + T, 0, 1]

>