# Suppose $g$ is a permutation of [1..n]

NrCyclesPerm := function(g, omega)


    local cs, i, t;

    if not IsSubset(omega,MovedPoints(g)) then return false; fi;

    cs := CycleStructurePerm(g);
    t := 0; # counts the number of cycles
    for i in [ 1..Length(cs) ] do
        if IsBound(cs[i]) then 
            t := t + cs[i]; 
        fi;
    od;

    # note that we have to add the points in 1-cycles
    return t + Length(omega) - Length(MovedPoints(g));
end;


NrColourings := function( grp, nr_col )

	local cc, c, a, t, omega;

    omega := MovedPoints(grp);
	# Compute the conjugacy classes of grp
	cc := ConjugacyClasses(grp);

    t := 0;
	# Loop over the conjugacy classes
	for c in cc do
		# choose a representative of class c
		a := Representative(c);
		t := t + nr_col^NrCyclesPerm(a,omega)/Size(Centralizer(grp,a));
	od;

	return t;

end;