anyone?

I can suggest you the following strightforward approach:

• derive or find the exact formula $$$S(n) = \Large \sum\limits_{d|n} \frac{n!}{d!\,(\frac{n}{d}!)^d}$$$
• represent large numbers as its prime factorisation; i've used

struct component
{
  short cnt[1229];
};

where 1229 is the number of prime numbers till 10000 ($$$N_{max} \times K_{max}$$$) and cnt is the power of corresponding prime number in order to transform all multiplications and divisions to additions and substractions

• use optimised bigint library for converting factoring numbers into digital form.