with(numtheory); pet_cycleind_cyclic := proc(n) local d; add(phi(d)*a[d]^(n/d), d in divisors(n))/n; end; pet_cycleind_dihedral := proc(n) local s; s := 1/2*pet_cycleind_cyclic(n); if(type(n, odd)) then s := s + 1/2*a[1]*a[2]^((n-1)/2); else s := s + 1/4*(a[1]^2*a[2]^((n-2)/2) + a[2]^(n/2)); fi; s; end; pet_cycleind_symm := proc(n) option remember; local l; if n=0 then return 1; fi; expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n)); end; # number of bracelets of B beads using at most C swappable colors bracelet_pg := proc(B, C) 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 B = 1 or C = 1 then return 1 fi; idx_slots := pet_cycleind_dihedral(B); idx_colors := pet_cycleind_symm(C); 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; bracelet_pg_OEIS := proc(B,C) if B=1 or C=1 then return 1 fi; bracelet_pg(B,C)-bracelet_pg(B,C-1); end;