Powerful Number 筛

本文主要讲解 Powerful Number 筛，并为之前的杜教筛进行一定的补充

$\text{powerful number}$ 定义：如果一个正整数 $x$ 所有质因子的幂次均大于等于二，即 $x=\prod p_i^{e_i}$ 且 $e_i\ge 2$ ，则称 $x$ 是一个 $\text{powerful number}$

一个重要的性质： $[1,n]$ 这个区间的 $\text{powerful number}$ 有 $O(\sqrt{n})$ 个
证明：
不难发现，所有的 $\text{powerful number}$ 都可以表示成 $a^2{b}^3$ 的形式。我们可以通过枚举 $a$ 来表示区间 $[1,n]$ 的 $\text{powerful number}$ 数量：
$$O(\sum _{a=1}^{\sqrt{n}}(\frac{n}{a^2}{)}^{\frac{1}{3}})=O({\int }_1^{\sqrt{n}}(\frac{n}{a^2}{)}^{\frac{1}{3}}da)=O(3n^{\frac{1}{3}}{a}^{\frac{1}{3}}{|}_1^{\sqrt{n}})=O(\sqrt{n})$$

Powerful Number 筛：
以 $P5325$ 为例，该题是求 $f(p^k)=p^k(p^k-1)$ 这个积性函数的前缀和
我们构造两个积性函数 $g(x),h(x)$ ，其满足 $g(p)=f(p)$ ， $f=g\ast h$
此时有 $f(p)=g(p)h(1)+h(p)g(1)$ ，而 $f(p)=g(p),h(1)=g(1)=1$ 所以 $h(p)=0$
而由于 $h(x)$ 是一个积性函数，所以其函数值不为 $0$ 的点自变量一定是 $\text{powerful number}$
那么题目所求可进行如下转化：
$$\sum _{i=1}^{n}f(i)=\sum _{ij\le n}g(i)h(j)=\sum _{i=1}^{n}h(i)\sum _{j=1}^{\lfloor \frac{n}{i}\rfloor }g(j)$$
Powerful Number 筛的瓶颈就在于求 ${\sum }_{j=1}^{\lfloor \frac{n}{i}\rfloor }g(j)$ ，对于一般的简单数论函数我们通常采取杜教筛来以 $O(n^{\frac{2}{3}})$ 来求 $g$ 的前缀和
对于本题 $f(p)=p(p-1)=p\phi (p)$ ，所以我们不妨令 $g(x)=x\phi (x)$
$g(x)$ 的前缀和我们可以用经典杜教筛求：
我们知道 $\phi \ast I=id$ ， 当 $C$ 是完全积性函数时有如下性质：
$(A\cdot C)\ast (B\cdot C)=(A\ast B)\cdot C$
证明：
$$\sum _{d|n}(A(d)C(d))(B(\frac{n}{d})C(\frac{n}{d}))=C(n)\sum _{d|n}A(d)B(\frac{n}{d})$$
所以我们可以通过给 $\phi \ast I=id$ 两边乘一个 $id$ 来构造杜教筛，即 $id\cdot \phi \ast id=i{d}^2$
这样我们就可以 $O(n^{\frac{2}{3}})$ 来求 ${\sum }_{i=1}^{n}i\phi (i)$ 了
代码如下：

ll get_sum(int x)
{
	if(x <= limit) return g[x];
	if(mp[x]) return mp[x];
	ll tmp = x%mod;
	ll res = x*(x+1)%mod*(2*x+1)%mod*inv6%mod;
	for(ll l = 2,r;l <= x;l = r+1)
	{
		r = x/(x/l);
		res = (res-get_sum(x/l)*(r-l+1)%mod*(l+r)%mod*inv2%mod)%mod;
	}
	return mp[x] = (res+mod)%mod;
}

接下来问题就只剩求 $h(x)$ ，因为 $\text{powerful number}$ 在 $[1,n]$ 只有 $O(\sqrt{n})$ 个，所以只需要求 $h$ 在 $\text{powerful number}$ 处的取值即可（其他地方的 $h(x)=0$）
这里有一种生成 $\text{powerful number}$ 的方法：
枚举 $p^e(e\ge 2,p\le \sqrt{n})$ 通过 $dfs$ 硬搞
对于 $h$ 我们考虑有迪利克雷卷积的定义式子，将其展开：
$$f(p^k)=\sum _{i+j=k}g(p^i)h(p^j)$$
所以移项有：
$$h(p^k)=f(p^k)-\sum _{i=0}^{k-1}g(p^{k-i})h(p^i)$$
然后利用积性函数的性质合并质数即可

具体细节见代码：

//#pragma GCC optimize(2)
//#pragma GCC optimize("Ofast","inline","-ffast-math")
//#pragma GCC target("avx,sse2,sse3,sse4,mmx")
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#include<set>
#include<map>
#include<stack>
#include<queue>
#include<unordered_map>
#define ll long long
#define inf 0x3f3f3f3f
#define int  ll
#define endl '\n'
#define IOS ios::sync_with_stdio(0); cin.tie(0); cout.tie(0)
using namespace std;
int read()
{
	int res = 0,flag = 1;
	char ch = getchar();
	while(ch<'0' || ch>'9')
	{
		if(ch == '-') flag = -1;
		ch = getchar();
	}
	while(ch>='0' && ch<='9')
	{
		res = (res<<3)+(res<<1)+(ch^48);//res*10+ch-'0';
		ch = getchar();
	}
	return res*flag;
}
const int maxn = 1e6+5;
const int mod = 1e9+7;
const double pi = acos(-1);
const double eps = 1e-8;
int n,cnt,limit,sq,pri[maxn],g[maxn],ans,h[maxn][35];
bool vis[maxn];
void init(int n)
{
	vis[1] = true; g[1] = 1;
	for(int i = 2;i <= n;i++)
	{
		if(!vis[i])
		{
			pri[++cnt] = i;
			g[i] = i-1;
		}
		for(int j = 1;j <= cnt && i*pri[j]<=n;j++)
		{
			vis[i*pri[j]] = true;
			if(i%pri[j] == 0) 
			{
				g[i*pri[j]] = g[i]*pri[j];
				break;
			}
			g[i*pri[j]] = g[i]*(pri[j]-1);
		}
		g[i] = (g[i-1]+g[i]*i)%mod;
	}
	for(sq = 1;pri[sq]*pri[sq] <= ::n;sq++);
	int gg[35];
	for(int i = 1;i <= sq;i++)
	{
		gg[0] = h[i][0] = 1;
		gg[1] = pri[i]*(pri[i]-1)%mod;
		int p2 = pri[i]*pri[i]%mod;
		for(int p = pri[i]*pri[i],e = 2;p <= ::n;p *= pri[i],e++)
		{
			gg[e] = gg[e-1]*p2%mod;
			h[i][e] = p%mod*(p%mod-1)%mod;
			for(int j = 0;j < e;j++) h[i][e] = (h[i][e]-h[i][j]*gg[e-j]%mod+mod)%mod;
		}
	}
} 
unordered_map<int,ll>mp;
ll qpow(ll a,ll b)
{
	ll res = 1;
	while(b)
	{
		if(b&1) res = res*a%mod;
		a = a*a%mod;
		b >>= 1;
	}
	return res;
}
const int inv2 = qpow(2,mod-2),inv6 = qpow(6,mod-2);
ll get_sum(int x)
{
	if(x <= limit) return g[x];
	if(mp[x]) return mp[x];
	ll tmp = x%mod;
	ll res = tmp*(tmp+1)%mod*(2*tmp+1)%mod*inv6%mod;
	for(ll l = 2,r;l <= x;l = r+1)
	{
		r = x/(x/l);
		res = (res-get_sum(x/l)*(r-l+1)%mod*((l+r)%mod)%mod*inv2%mod)%mod; //注意 l+r 取模 
	}
	return mp[x] = (res+mod)%mod;
}
void dfs(int cur,int val,int st)
{
	ans = (ans+val*get_sum(n/cur))%mod;
	for(int i = st;i <= sq;i++)
	{
		if(pri[i]*pri[i] > n/cur) return ;
		for(ll p = pri[i]*pri[i],e = 2;p*cur <= n;p *= pri[i],e++)
			dfs(cur*p,val*h[i][e]%mod,i+1);
	}
}
signed main()
{
	n = read();
	limit = max((int)pow(n,0.6),100LL);
	init(limit);
	dfs(1,1,1);
	cout<<ans<<endl;
	return 0;
}