// #pragma GCC optimize("O3")
// #pragma GCC optimize("Ofast")
// #pragma GCC optimize("unroll-loops")
// #pragma GCC target("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx,avx2,tune=native")
#include "pilgrimage.h"
#include <bits/stdc++.h>
using namespace std;

#define PB emplace_back
#define int long long
#define ll long long
#define vi vector<int>
#define siz(a) ((int) ((a).size()))
#define rep(i, a, b) for (int i = (a); i <= (b); ++i)
#define per(i, a, b) for (int i = (a); i >= (b); --i)

const int N = 1e6;
struct fastmod {
	typedef unsigned long long u64;
	typedef __uint128_t u128;
	int m;
	u64 b;
	fastmod(int m) : m(m), b(((u128)1 << 64) / m) {}
	int reduce(u64 a) {
		u64 q = ((u128)a * b) >> 64;
		int r = a - q * m;
		return r < m ? r : r - m;
	}
} z(2);
struct pii {int x, y;};
namespace solve {
	ll n, m, p;
	inline void Exgcd(ll a, ll b, ll &x, ll &y) {
		if (!b) x = 1, y = 0;
		else Exgcd(b, a % b, y, x), y -= a / b * x;
	}
	inline ll fpm(ll x, ll power, ll Mod) {
		ll res = 1;
		for (; power; power >>= 1, (x *= x) %= Mod)
			if (power & 1) (res *= x) %= Mod;
		return res;
	}
	inline ll fac(ll n, ll pi, ll pk) {
		if (!n) return 1;
		ll res = 1;
		rep (i, 2, pk) if (i % pi) (res *= i) %= pk;
		res = fpm(res, n / pk, pk);
		rep (i, 2, n % pk) if (i % pi) (res *= i) %= pk;
		return res * fac(n / pi, pi, pk) % pk;
	}
	inline ll Inv(ll n, ll Mod) {
		ll x, y;
		Exgcd(n, Mod, x, y);
		return (x % Mod + Mod) % Mod;
	}
	inline ll CRT(ll b, ll Mod) {
		return b * Inv(p / Mod, Mod) % p * (p / Mod) % p;
	}
	inline ll factor(ll x, ll Mod) {
		return x ? factor(x / Mod, Mod) + (x / Mod) : 0;
	}
	inline ll Comb(ll n, ll m, ll pi, ll pk) {
		ll k = factor(n, pi) - factor(m, pi) - factor(n - m, pi);
		if (!fpm(pi, k, pk)) return 0;
		return fac(n, pi, pk) * Inv(fac(m, pi, pk), pk) % pk * Inv(fac(n - m, pi, pk), pk) % pk * fpm(pi, k, pk) % pk;
	}
	vector<pii> s;
	inline ll ExLucas(ll n, ll m) {
		ll res = 0, tmp = p;
		if(!siz(s)) {
			rep(i, 2, sqrt(p + .5)) if (!(tmp % i)) {
				ll pk = 1;
				while (!(tmp % i)) pk *= i, tmp /= i;
				res = z.reduce(res + CRT(Comb(n, m, i, pk), pk));
				s.PB((pii) {i, pk});
			}
			if (tmp > 1) res = z.reduce(res + CRT(Comb(n, m, tmp, tmp), tmp)), s.PB((pii) {tmp, tmp});
		}
		else {
			for(pii x : s) res = z.reduce(res + CRT(Comb(n, m, x.x, x.y), x.y));
		}

		return res;
	}
	int ANS (int _n, int _m, int _p) {
		n = _n, m = _m, p = _p;
		return ExLucas(n, m);
	}
}

int a, b, i2, mod;
int qp(int n, int m = mod - 2) {
	int res = 1;
	for(; m; m >>= 1) {
		if(m & 1) res = z.reduce(res * n);
		n = z.reduce(n * n);
	}
	return res;
}
int phi(int n) {
	int res = n;
	for(int i = 2; i * i <= n; ++i) {
		if(n % i == 0) {
			res = res / i * (i - 1);
			for(; n % i == 0; n /= i);
		}
	}
	if(n > 1) res = res / n * (n - 1);
	return res;
}

#undef int
void init (int o, int p) {
	z = fastmod(p);
	mod = p;
	a = phi(p) - 1;
	i2 = qp(2, a);
}
int ask(ll n) {
	// cerr << solve::ANS(2 * n - 1, n, mod) << endl;
	return solve::ANS(2 * n - 1, n, mod) * (n % mod) * 2 % mod * qp(i2, 2 * n) % mod;
}

// signed main() {
// 	// freopen("in1.in", "r", stdin);
// 	// freopen("out.out", "w", stdout);
// 	// cout << calc(666, 64, 2).val;
// 	// return 0;
// 	printf("%lld", ans);
// 	return 0;
// }