參考部落格：

https://www.cnblogs.com/zzqsblog/p/9904271.html

powerful number：
每個質因子指數\(\ge 2\)的數。

可以被\(a^2 b^3(\mu(b) \not= 0)\)唯一表示。

也可以被\(ab^2(a|b)\)唯一表示。

題外（powerful number計數）

https://loj.ac/article/589

\(n\)以內的powerful number的數量大概是\(2\sqrt n\)

利用powerful number篩積性函式字首和

https://gmoj.net/senior/#main/show/6785

例子：
設積性函式\(f(x)=2^{x的指數和}\)

設\(d(x)\)表示\(x\)的約數個數。

我們觀察\(g=f/d\)是什麼：

因為都是積性函式，所以可以只考慮\(p^k\)的時候：
以下為\(k=0,1,2,3,…\)的情況：
\(f=1,2,4,8,…\)
\(d=1,2,3,4,…\)
\(g=1,0,1,2,…\)

不難發現\(g(p^1)=0\)，也就是\(g\)只在powerful number處有值，考慮反推\(f\)。

\(f=g*d\)
\(\sum_{i=1}^n f(i)=\sum_{x~is~a~powerful~number} g(x)*\sum_{j=1}^{n/x} d(j)\)

可以線性篩\(1-\sqrt n\)的約數個數，這樣\(n/x\le sqrt(n)\)時直接查詢，否則利用\(\sum_{i=1}^n d(i)=\sum_{i=1}^n n/i\)分塊解決。

通過積分可分析出複雜度為\(O(\sqrt n log n)\)。

事實上還可以解決下面積性函式：

1. \(f(p)=1\)的，使\(f\)除以\(id0\)函式即可。
2. \(f(p)=p^k\)的，使\(f\)除以\(idk\)函式即可。
3. \(f(p)=p-1\)的，使\(f\)除以\(\phi\)函式即可。
4. \(f(p)=2\)的，使\(f\)除以\(id0*id0\)（約數個數）函式即可。
5. \(f(p)=p+1\)
   的，使\(f\)除以\(id0*id1\)（約數和）函式即可。

不難發現只要除的函式是個字首和好算的就行了。

其實運用還是很侷限的，

參考程式碼：

#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, _b = y; i <= _b; i ++)
#define ff(i, x, y) for(int i = x, _b = y; i <  _b; i ++)
#define fd(i, x, y) for(int i = x, _b = y; i >= _b; i --)
#define ll long long
#define pp printf
#define hh pp("\n")
using namespace std;

ll n, mo, sq;

ll ksm(ll x, ll y) {
	ll s = 1;
	for(; y; y /= 2, x = x * x % mo)
		if(y & 1) s = s * x % mo;
	return s;
}

ll a2[125], xs[125];

void build(int n) {
	a2[0] = 1; fo(i ,1, n) a2[i] = a2[i - 1] * 2 % mo;
	fo(i, 0, n) {
		xs[i] = a2[i];
		fo(j, 0, i - 1) xs[i] = (xs[i] - xs[j] * (i - j + 1)) % mo;
	}
}

const int N = 1e7 + 5;

int bz[N], p[N], p0, e[N];
ll cd[N];

void sieve(int n) {
	cd[1] = 1;
	fo(i, 2, n) {
		if(!bz[i]) {
			p[++ p0] = i;
			e[i] = 1;
			cd[i] = 2;
		}
		for(int j = 1; i * p[j] <= n; j ++) {
			int k = i * p[j]; bz[k] = 1;
			if(i % p[j] == 0) {
				e[k] = e[i] + 1;
				cd[k] = cd[i] / (e[i] + 1) * (e[k] + 1);
				break;
			}
			e[k] = 1;
			cd[k] = cd[i] * 2;
		}
	}
	fo(i, 1, n) cd[i] = (cd[i - 1] + cd[i]) % mo;
}

ll ans = 0;

int cnt = 0;

void tj(ll xs, ll x) {
	ll m = n / x;
	if(m <= sq) {
		ans = (ans + cd[m] * xs) % mo;
	} else {
		ll sum = 0;
		ll c = sqrt(m);
		for(ll i = 1; i <= c; i ++) {
			sum += m / i;
		}
		sum = sum * 2 - c * c;
		ans = (ans + sum % mo * xs) % mo;
	}
}

void dg(int i, ll s, ll s2) {
	tj(s2, s);
	if(i > p0) return;
	if(n / p[i] / p[i] < s) return;
	for(int j = i; j <= p0 && (ll) s * p[j] * p[j] <= n; j ++) {
		ll p1 = p[j];
		for(int e = 2; n / p[j] / p1 >= s; e ++) {
			ll p2 = p1 * p[j];
			dg(j + 1, s * p2, s2 * xs[e]);
			p1 = p2;
		}
	}
}

int main() {
	freopen("remapping.in", "r", stdin);
	freopen("remapping.out", "w", stdout);
	scanf("%lld %lld", &n, &mo)	;
	sq = sqrt(n);
	build(120);
	sieve(sq);
	dg(1, 1, 1);
	pp("%lld\n", ans);
}