利用powerful number筛积性函数前缀和

参考博客：

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的指数和}\)
求\(\sum_{i=1}^n f(i)(n=10^{14})\)

设\(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);
}