cometoj 1098 神奇函数【数论+powerful number筛】

cometoj 1098 【数论+powerful number筛】

”即便它再像一个能筛的函数，也要先推它的定义式“。——未来哲学家yxy

若 $x=\prod p_i^{e_i}$
$$\begin{array}{rl}f(x)& =\prod p_i^{\lfloor e^i/2\rfloor }\\ & =\frac{x}{\prod p_i^{\lceil e^i/2\rceil }}\\ & =\sum _{j=1}^x[\prod _{i=1}^k{p}_i^{\lceil e^i/2\rceil }|j]\\ & =\sum _{j=1}^x[x|j^2]\end{array}$$
第二部到第三步用了的除法的定义 $\frac{a}{b}={\sum }_{i=1}^a[b|i]$ ，意思是 $a$ 里面有多少个 $b$ 的倍数。

最后一步意思是当 $a|b$ ，当且仅当 $a^2|b^2$ （就**离谱。

考虑求前缀和，我们知道 $f(x)$ 是积性函数，且 $f(p)=1,p\in Prime$ 。

powerful number筛，结束。

#include <bits/stdc++.h>
#define N 3200006
using namespace std;
typedef unsigned long long ll;
ll n,sqr,pri[N],tot,ans;
bool tag[N];
void init(){
	sqr=sqrt(1e13); 
	for(ll i=2;i<=sqr;i++){
		if(!tag[i]) pri[++tot]=i;
		for(int j=1;pri[j]*i<=sqr;j++){
			tag[pri[j]*i]=1;
			if(i%pri[j]==0)break;
		}
	}
}
void solve(int pos,ll now,ll lst){
	ll lim=n/now;
	if(pos==tot+1||now*pri[pos]>n||now*pri[pos]*pri[pos]>n){
        ans+=lst*(n/now); return;
    }
	ll p=pri[pos],nxt=1,pp=1; bool h=1;
	solve(pos+1,now,lst);
	while(1){
		if(nxt*p>lim)break;nxt*=p;
		if(nxt*p>lim)break;nxt*=p;
		if(h)lst*=p-1,h=0; else lst*=p;
		solve(pos+1,nxt*now,lst);
	}
}
int main(){
	int T; cin>>T;
	init(); 
	while(T--){
		cin>>n; 
		solve(1,1,1);
		cout<<ans<<'\n'; ans=0;
	}
}