M := proc(n,b) local dv,s,d;
 dv := numtheory[divisors](n);
 s := 0;
 for d in dv do
  s := s+numtheory[mobius](iquo(n,d))*((b+1)^d-b^d);
 od;
 s/n;
end:

M(3,3);

### For use by authors only
### Experimental editor for the Maple Form Interface
###