(* pu3.m:  S. A. Fulling, Nov. 1995 (first version summer 1994) *)

Off[General::spell1]
Share[]
orderRules = {
   x1^(n1_)*x2^(n2_) :> x1^n2*x2^n1 /; n2 > n1,
   x1^(n1_)*x3^(n3_) :> x1^n3*x3^n1 /; n3 > n1,
   x2^(n2_)*x3^(n3_) :> x2^n3*x3^n2 /; n3 > n2
   }

qmax = 8
max = 2(qmax+1)

link12 = Sum[t^(2j) c[j] x1^j x2^j, {j, 0, qmax}] +O[t]^max
link13 = link12 /. x2 -> x3
link23 = link13 /. x1 -> x2
point1 = Sum[t^(2j) c[j] x1^(2j), {j, 0, qmax}] + O[t]^max
point2 = point1 /. x1 -> x2
point3 = point1 /. x1 -> x3
pointpair = Sum[t^(4j) c[j]^2 x1^(2j) x2^(2j), 
                                {j, 0, Floor[qmax/2]}] + O[t]^max
angle = Sum[t^(4j) c[j]^2 x1^j x2^j x3^(2j), 
                                {j, 0, Floor[qmax/2]}] + O[t]^max
pointtriple = Sum[t^(6j) c[j]^3 x1^(2j) x2^(2j) x3^(2j), 
                                {j, 0, Floor[qmax/3]}] + O[t]^max
linktriple = pointtriple

identity = link12*link13*link23*point1*point2*point3
transpos = link12*angle*pointpair*point3
cyclics = pointtriple*linktriple

pu3series = (1/6)(identity + 3*transpos + 2*cyclics)
pu3temp = (Expand[Normal[pu3series] * x1^2 x2^2 x3^2] //. 
                                 orderRules) / (x1^2 x2^2 x3^2)
pu3 = Collect[pu3temp, Prepend[Table[c[qmax-i], {i,0,qmax}], t]]

pu3point = pu3 //. c[n_] -> 1
pu3line = pu3 //. {x1->1, x2->1, x3->1}
pu3total = pu3line //. c[n_] -> 1

Save["pu3.out", pu3, pu3point, pu3line, pu3total]
Quit[]