with(numtheory); pet_cycleind_square := proc(n) option remember; local cind; if type(n, even) then cind := 1/8*(a[1]^(n^2) # identity +2*a[2]^(n^2/2) # horiz. / vert. flips +2*a[1]^n*a[2]^((n^2-n)/2) # diagonal flips +a[2]^(n^2/2) # 180 deg. rotation +2*a[4]^(n^2/4)); # 90 / 270 deg. rotation else cind := 1/8*(a[1]^(n^2) # identity +2*a[1]^n*a[2]^((n^2-n)/2) # horiz, / vert. flips +2*a[1]^n*a[2]^((n^2-n)/2) # diagonal flips +a[1]*a[2]^((n^2-1)/2) # 180 deg. rotation +2*a[1]*a[4]^((n^2-1)/4)); # 90 / 270 deg. rotation fi; return cind; end; pet_cycleind_symm := proc(n) option remember; if n=0 then return 1; fi; expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n)); end; square_pg := proc(side, colors) option remember; local idx_slots, idx_colors, res, term_a, term_b, v_a, v_b, inst_a, inst_b, len_a, len_b, p, q; if side = 1 or colors = 1 then return 1 fi; idx_slots := pet_cycleind_square(side); idx_colors := pet_cycleind_symm(colors); res := 0; for term_a in idx_slots do for term_b in idx_colors do p := 1; for v_a in indets(term_a) do len_a := op(1, v_a); inst_a := degree(term_a, v_a); q := 0; for v_b in indets(term_b) do len_b := op(1, v_b); inst_b := degree(term_b, v_b); if len_a mod len_b = 0 then q := q + len_b*inst_b; fi; od; p := p*q^inst_a; od; res := res + lcoeff(term_a)*lcoeff(term_b)*p; od; od; res; end;