When I was making a problem, I needed to calculate $$$C_n^r$$$ modulo a non-prime big number. So, I found a simpler method I haven't seen shared, and I'm excited to explain it.

Prerequisites:

• Modular arithmetic

• Factorization

• Binary exponentiation

• Knowing $$$\binom{n}{r}$$$ and Euler's totient function

The problem is the modular inverse. We can get it for a number $$$x$$$ if and only if $$$x$$$ is coprime with the modulus. But in a lot of problems, we just need to undo a previous operation,

$$$C_5^3 = \frac{5!}{3! \cdot 2!}$$$. We delete one $$$3$$$ and two $$$2$$$s from $$$5!$$$.

So why have we never thought to store it in a frequency array?

My idea is that when multiplying, we use the modular arithmetic formulas, and when dividing, we will count the frequency of each prime number in the modulus and delete it from $$$n!$$$, So we need to first calculate $$$n!$$$ without the primes in the modulus and count the frequency for each prime. Then, when you need to know the answer, you can iterate over the frequency array :)

But what if I divided $$$\frac{1}{2}$$$? It will not even happen, so what is the meaning of $$$\frac{1}{2}$$$?

I think it looks like imaginary numbers. $$$\sqrt{−1}^2 = −1$$$. It's easy to understand :)

We will find the prime factors of a number with any method (normal factorization is good enough). The modular inverse for a coprime number is $$$x^{phi(mod) - 1}$$$ (it's not working in power towers, but here I think it's enough here), $$$phi(x)$$$ is a multiplicative function you can easily calculate in $$$O(\sqrt{n})$$$.

So we can imagine a number, let's call him imaginary-modular-integer (imint), and make the necessary operations only (I mean * and /, no need for *= and /=).

Note: if you need to divide after adding, my trick will never work unless you make it better :) (tell me if you did).

You can see the code here:

struct imint {
    std::vector<__uint64_t>_primes, _power;
    __uint64_t _modulo, _phi;
    __uint128_t _value;
    void build(__uint64_t _target) {
        _modulo = _target, _phi = _modulo, _value = 1;
        for (__uint64_t i = 2; i * i <= _modulo; ++i)
        if (_modulo % i == 0) {
            while (_modulo % i == 0) _modulo /= i;
            _phi -= _phi / i;
        }
        if (_modulo > 1) _phi -= _phi / _modulo;
        _modulo = _target;
        for (__uint64_t i = 2; i * i <= _target; ++i)
            if (_target % i == 0) {
                __uint64_t j = _target / i;
                _primes.push_back(i);
                if (i * i != _target)
                    _primes.push_back(j);
                while (_target % i == 0) _target /= i;
                while (_target % j == 0) _target /= j;
            }
            if (_target > 1) _primes.push_back(_target);
            _power.assign(_primes.size(), 0);
    }
    friend imint operator*(__uint64_t a, imint b) {
        for (size_t i = 0; i < b._power.size(); i++)
            while (a % b._primes[i] == 0) {
                b._power[i]++;
                a /= b._primes[i];
            }
        (b._value *= a) %= b._modulo;
        return b;
    }
    friend imint operator/(imint a, imint b) {
        assert(a._modulo == b._modulo);
        for (size_t i = 0; i < a._power.size(); ++i)
            a._power[i] -= b._power[i];
        __uint64_t _exponentiation = a._phi - 1;
        while (_exponentiation) {
            if (_exponentiation & 1)
                (a._value *= b._value) %= a._modulo;
            (b._value *= b._value) %= a._modulo;
            _exponentiation >>= 1;
        }
        return a;
    }
    __uint64_t real() {
        __uint128_t _real = _value;
        for (size_t i = 0; i < _power.size(); ++i) {
            __uint32_t _exponentiation = _power[i];
            __uint128_t _base = _primes[i];
            while (_exponentiation) {
                if (_exponentiation & 1)
                    (_real *= _base) %= _modulo;
                (_base *= _base) %= _modulo;
                _exponentiation >>= 1;
            }
        }
        return _real;
    }
};

After using this code as a template (please understand it first, don't just copy it!), you can write this:

dp[0].build(m);
for (int i = 1; i <= n; i++) dp[i] = i * dp[i - 1];

It will add $$$\lg(\lg(n))$$$ to the time complexity :(

When you need to calculate the answer, you can use $$$real()$$$ function.

Important Note: Do not write

dp[n] / (dp[r] * dp[n - r]);

. Instead, you should use:

(dp[n] / dp[r] / dp[n - r]).real();

It's not hard and pretty cool. If you know stars and bars and sieve, you can try to solve the easy version of my problem: Invitation

These are hardest two problems from Today's contest in Asyut Hack Club :)